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Abstract 

The distribution of organisms in space has important consequences for the function 
and structure of ecological systems. Such distributions are often referred to as patchy, 
and a patch-based approach to modeling ecosystem dynamics has become a major 
research focus. These models have been used to explore a wide range of questions 
concerning population, metapopulation, community, and landscape ecology, in both 
terrestrial and aquatic systems. 

In this dissertation I develop and analyze a series of spatial models to study the dy¬ 
namics of metapopulations and marine benthic communities in patchy environments. 
All the models have the form of a discrete-time Markov chain, and assume that the 
landscape is composed of discrete patches, each of which is in one of a number of 
possible states. The state of a patch is determined by the presence of an individual 
of a given species, a local population, or a group of species, depending on the spatial 
scale of the model. 

The research is organized into two main parts as follows. In the first part, I 
present an analysis of the effects of habitat destruction on metapopulation persistence. 
Theoretical studies have already shown that a metapopulation goes extinct when the 
fraction of suitable patches in the landscape falls below a critical threshold (the so 
called extinction threshold). This result has become a paradigm in conservation 
biology and several models have been developed to calculate extinction thresholds 
for endangered species. These models, however, generally do not take into account 
the spatial arrangement of habitat destruction, or the actual size of the landscape. 
To investigate how the spatial structure of habitat destruction affects persistence, 
I compare the behavior of two models: a spatially implicit patch-occupancy model 
(which recreates the extinction patterns found in other models) and a spatially explicit 
cellular automaton (CA) model. In the CA, I use fractal arrangements of suitable 
and unsuitable patches to simulate habitat destruction and show that the extinction 
threshold depends on the fractal dimension of the landscape. To investigate how 
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habitat destruction affects persistence in finite landscapes , I develop and analyze 
a chain-binomial metapopulation (CBM) model. This model predicts the expected 
extinction time of a metapopulation as a function of the number of patches in the 
landscape and the number of those patches that are suitable for the population. 
The CBM model shows that the expected time to extinction decreases greater than 
exponentially as suitable patches are destroyed. I also describe a statistical method 
for estimating parameters for the CBM model in order to evaluate metapopulation 
viability in real landscapes. 

In the second part, I develop and analyze a series of Markov chain models for 
a rocky subtidal community in the Gulf of Maine. Data for the model comes from 
ten permanent quadrats (located on Ammen Rock Pinnacle at 30 meters depth) 
monitored over an 8-year period (1986-1994). I first parameterize a linear (ho¬ 
mogenous) Markov chain model from the data set and analyze it using an array of 
novel techniques, including a compression algorithm to classify species into functional 
groups, a set of measures from stochastic process theory to characterize successional 
patterns, sensitivity analyses to predict how changes in various ecological processes 
effect community composition, and a method for simulating species removal to iden¬ 
tify keystone species. I then explore the effects of time and space on successional 
patterns using log-linear analysis, and show that transition probabilities vary sig¬ 
nificantly across small spatial scales and over yearly time intervals. I examine the 
implications of these findings for predicting equilibrium species abundances and for 
characterizing the transient dynamics of the community. Finally, I develop a non¬ 
linear Markov chain for the rocky subtidal community. The model is parameterized 
using ma ximum likelihood methods to estimate density-dependent transition prob¬ 
abilities. I analyze the best fitting models to study the effects of nonlinear species 
interactions on community dynamics, and to identify multiple stable states in the 
subtidal system. 


Dissertation Advisor: Hal Caswell 

Title: Senior Scientist, Woods Hole Oceanographic Institution 
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Chapter 1 


Introduction 


Here is the world, sound as a nut, perfect, not the smallest piece of chaos 
left, never a stitch nor an end, not a mark of haste, or botching, or second 
thought; but the theory of the world is a thing of shreds and patches. 

-Ralph Waldo Emerson 

A characteristic feature of the spatial distribution of most species, across a range 
of spatial scales, is patchiness (Taylor 1961, Taylor and Taylor 1969; Hanski 1994; 
Wiens 1997). Patchiness refers to the non-homogenous distribution of organisms (and 
their resources) in space and time (Hanski 1999). At smaller spatial scales, patchiness 
is generally a consequence of neighborhood effects; i.e., organisms are more likely to 
interact with their neighbors than with more distant organisms. This is especially 
true for benthic invertebrates such as corals, sponges, bryozoans, and other sessile 
marine organisms (Tanner et al. 1994, Sebens 1986, Witman and Dayton 2000). At 
larger spatial scales, patchiness arises due to dispersal and recruitment patterns (e.g. 
Horn and MacArthur 1972, Levin 1974, Levin et at. 1984, Paine 1984, Cohen and 
Levin 1991, Caswell and Cohen 1991a,b, Hanski 1999), physical disturbance (Dayton 
et al 1970, Witman 1987, Witman and Dayton 2000), and variations in the quality 
of the environment (Kareiva 1990; Levin 1992; Tilman and Kareiva 1997). Because 
the absence of a species from a locality may reflect dispersal limitations, unsuitable 
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environmental conditions, or competitive exclusion through local biotic interactions, 
patchy distributions are an emergent property of ecological processes acting over a 
wide range of spatial scales (e.g., Horn 1971, Hastings 1980, Connell 1985, Gaines 
and Roughgarden 1985, Menge and Sutherland 1987, Grosberg and Levitan 1992, 
Connolly and Roughgarden 1998). 

The importance of spatial patchiness is increasingly recognized as having impor¬ 
tant consequences for the function and structure of terrestrial and aquatic systems 
(Levin and Paine 1974; Steel 1978; Pickett and White 1985; Caswell and Cohen 
1991a,b; Weins et. al. 1993; Wu and Levin 1994; Wu and Loucks 1995; Hanski 
1999). A patch is generally defined in terms of a spatial unit within a landscape that 
is different in some physical or biological aspect from its immediate surroundings 
(Kotliar and Weins 1990). Definitions of a patch include a bounded region within a 
homogenous background consisting of either single or multiple biological components 
(Levins and Paine 1974); a spatial location within an environment where the abun¬ 
dance of a resource, a population, or a community of organisms is high (Roughgarden 
1977); a location containing an aggregated collection of prey species within which a 
predator forages (Stephens and Kerbs 1986); territorial sites of individuals (Lande 
1987); discrete regions in a landscape containing local populations connected by dis¬ 
persal (Hanski 1994); and any division or heterogeneity in the abundance of resources 
(Antolin and Addicott 1991). Thus what is considered a patch depends both on the 
spatial and temporal scales of interest and the fundamental units of the system being 
studied (e.g., vertical nutrient distributions at the sediment water interface, territo¬ 
ries of spotted owls in fragmented forests, aggregates of benthic organisms, the global 
distribution of plankton). Spatial patchiness at any scale, however, can be defined 
in terms of both the composition of patches (i.e. types of patches and their relative 
abundance) and the spatial distribution of patches (i.e. patch size, shape, and their 
location in space). 

A patch-based approach to studying ecosystem dynamics has become a major 
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research focus, and several theoretical methods have been developed to model patch 
dynamics (Levins 1969; Lande 1989; Caswell and Cohen 1991a,b, 1995; Levin 1992; 
Weins et. al. 1993; Han ski 1994; Wu and Loucks 1995; Barradas et al. 1996; 
Caswell and Etter 1999). They include metapopulation models (Levin 1969, Nee 
and May 1992; Hanski 1999), patch-occupancy models (Caswell and Cohen 1991a,b, 
1995; Barradas et al. 1996), reaction-diffusion networks (Karlin and McGregor 1972; 
Levin 1974; Hastings 1983), Markov chains (Horn 1975; Usher 1979; Tanner et al. 
1994,1996), coupled map lattices (Comins et al. 1992; Sole and Vails 1991; Bevers 
and Flather 1999), interacting particle systems (Maekay and Jan 1994; Sutherland 
and Jacobs 1994; Durrett and Levin 1994), and cellular automata (Sivlertown et al. 
1992; Caswell and Etter 1993, 1999; Molofsky 1994; Dythan 1995; Bascompte and 
Sole 1996). These modeling approaches have been developed to answer a wide range 
of questions concerning population ecology, metapopulation dynamics, community 
ecology, biogeography, and landscape ecology (Wu and Loucks 1995). 

In this dissertation I develop a series of models to study the effects of habitat 
destruction on metapopulation persistence, and to analyze the dynamics of marine 
benthic communities. The general methods I use to model these systems all have 
the form of a discrete-time Markov chain. The models assume that the landscape 
is composed of discrete patches, each of which is in one of a number of possible 
states. The state of a patch is determined by the presence of an individual of a given 
species, a local population, or a group of species, depending on the spatial scale of 
the model. Before describing the specific research presented in the coming chapters, I 
briefly review the Markov chain approach used to model biological systems in patchy 
environments. 
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1.1 Characterizing patch dynamics: A discrete¬ 
time Markov chain approach 

Consider a landscape composed of patches, each of which is in one of N possible states. 
Let Xi(t) be the proportion of patches in state i at time t. The temporal dynamics of 
the patch states in the landscape can be modeled as a system of difference equations 

xi(t + 1) = fi(x 1 (t),x 2 {t),...,x N {t)) 
x 2 {t + l) = f 2 (xi(t),x 2 {t),...,x N (t)) 

i (1.1) 

x N (t + 1) = f N (xi(t),x 2 (t),...,x N (t)) 

where the functions /* specify the change in the proportion of patches in state i in the 
interval (t,t+ 1). If we define a vector x(t) whose ith element is £j(t), then equation 
1.1 can be written in matrix form 

x(t + 1) = Ax(t), (1.2) 

where A is a Markov chain transition matrix whose (i,j) element (a^) gives the prob¬ 
ability that a patch in state j changes to state i in one time interval. The matrix A 
is nonnegative (all ay > 0) and has the property that each column sums to 1 (i.e. A 
is column-stochastic). 

The set of possible state transitions for a given biological system can be depicted 
graphically using a transition diagram (Caswell and Cohen 1991a,b). Figure 1.1 shows 
an example of a transition diagram for a landscape in which there are four possible 
patch states. The arrows show the possible transitions among the patch states, while 
the expressions above the arrows are the probabilities that the transitions occur in the 
time interval (t, t + 1). Transition probabilities can be constant (e.g., S), a function 
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of time (e.g., 7 (t)), and/or dependent on the frequency of patches in a particular 
state (e.g., C'i(x 2 )). To construct a transition matrix from the diagram simply set the 
elements equal to the expression on the arrow going from X 3 to Xi. The transition 
matrix corresponding to Figure 1.1 is 

1 — Ci (x 2 ) 5 0 a 

Cx{x x ) {1 — C 2 (x 3 )}{1 — <5} 0 0 

0 C 2 (x 3 ){ 1 - <5} 1 - 7 (*) 0 

0 0 7(f) 1 — a 

1.2 Thesis structure 

The research presented in the following chapters focuses on the use of Markov chain 
models to investigate patch dynamics at the level of the metapopulation and at the 
scale of the community. While these models are applicable to many groups of organ¬ 
isms, they are primarily motivated by marine benthic systems. 

The models of metapopulation dynamics focus on the effects of habitat destruction 
on persistence. Theoretical models have shown that a metapopulation cannot persist 
when the fraction of suitable patches in the landscape falls below a critical threshold 
(Lande 1987; Nee and May 1992; Lamerson et al. 1994; Kareiva and Wennergren 
1995; Noon and McKelvey 1996). Consequently, identifying extinction thresholds for 
endangered species has become an important paradigm in conservation biology. These 
studies, however, have mostly ignored two fundamental factors affecting metapopu¬ 
lation persistence: ( 1 ). the spatial arrangement of habitat destruction, and ( 2 ). the 
actual size of the patch network. In the first part of the thesis (Chapters 2 and 
3 ) I develop and analyze a set of models that examine how changes in the habitat 
destruction pattern and changes in the size of the patch network affect predictions 
about metapopulation persistence and extinction thresholds in fragmented habitats. 
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1 —C,(x 2 ) 


(l-C 2 (X 3 )Xl-5) 



l-o 


i-Y(t) 


Figure 1.1: Transition diagram for a four state Markov chain model. Each circle 
represents a possible state of a patch, and in this example the states and are identified 
as 1 2, 3, and 4. The state of a patch is defined by the presence and/or absence of a 
particular species or species group. For instance, if the diagram were for a predator- 
prey model then 1 = empty patches, 2 = prey only, 3 = prey and predators, and 4 = 
predators only. The arrows connecting the circles represent possible state transitions, 
which the coefficients above the arrows are the transition probabilities during the 
time period (t, t + 1). 
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In second part of the thesis (Chapters 4, 5 and 6), I focus on the construction and 
analysis of a series of Markov chain models for a rocky subtidal community in the Gulf 
of Maine. The purpose of this research is four-fold. First, it extends the use of Markov 
chain models as a tool for marine community ecology. Second, it characterizes how 
temporal and spatial variation in successional processes affects species abundances 
and co mmuni ty dynamics. Third, it provides the first evaluation of nonlinear Markov 
models developed from empirical data. Finally, it increases our understanding of the 
ecology of sessile invertebrate communities in the rocky subtidal zone, a vast habitat 
that faces an increasing pace of anthropomorphic disturbance. 

1.2.1 Chapter details 

The chapters in this thesis are organized as follows: 

• In chapter 2 I explore the effect of habitat fragmentation when the fragmenta¬ 
tion follows a fractal pattern. The goals are to determine the habitat fragmen¬ 
tation pattern affects metapopulation abundance and the amount of habitat 
destruction the metapopulation can tolerate before it goes extinct. To study 
these effects, I compare the behavior of two models: a spatially implicit patch- 
occupancy (PO) model and a spatially explicit cellular automaton (CA) model. 

The PO model is a nonlinear Markov chain that describes patch dynamics 
within a fragmented habitat, in which patches are either suitable (i.e. can 
support local population growth) or unsuitable. Suitable patches are either 
occupied or unoccupied, and change state depending on rates of colonization 
and local extinction. The PO model recreates the extinction patterns found in 
other metapopulation models (e.g., Lande 1987, Hanski 1999). The advantage 
of the PO model is that it can be directly translated into a stochastic CA model 
(Caswell and Etter 1993, 1999), in which the state and the location of patches 
are followed explicitly through time. By using fractal arrangements of suitable 
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and unsuitable patches to simulate habitat fragmentation, I show that landscape 
structure plays a vital role in determining the effects of habitat destruction on 
persistence. 

• In chapter 3 I present a stochastic model for metapopulations in landscapes 
with a finite but arbitrary number of patches. The model, similar in form to 
chain-binomial epidemic models, is an absorbing Markov chain that describes 
the changes in the number of occupied patches as a sequence of binomial prob¬ 
abilities. It predicts the quasi-equilibrium distribution of occupied patches, the 
expected extinction time, and the probability of persistence to a given age as 
a function of the number of patches in the landscape and the number of those 
patches that are suitable for the population. In this chapter I also describe a sta¬ 
tistical method for estimating model parameters from time series data in order 
to evaluate metapopulation viability in real landscapes. An example is pre¬ 
sented using published data on the Glanville fritillary butterfly Meltiaea cinxia 
and its specialist parasitoid Cotesia melitaearum in which the expected extinc¬ 
tion time of M. cinxia is calculated as a function of the frequency of parasite 
outbreaks. 

• In chapter 4 I describe the development of a homogenous Markov chain model 
for a rocky subtidal community in the Gulf of Maine. Data for the model comes 
from ten permanent quadrats (located on Ammen Rock Pinnacle at 30 meters 
depth) monitored over an 8-year period (1986-1994). I analyze this model using 
an array of novel techniques, including a compression algorithm for transition 
matrices to classify species into functional groups, a set of tools from stochastic 
process theory to characterize successional patterns, and sensitivity analyses to 
quantify the importance of various ecological processes for maintaining species 
diversity. 

I also present methods for quantifying the effects of species removal on species 
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diversity and community resilience. These methods are Markov chain analogs 
to species removal experiments conducted in natural systems. They provide a 
theoretical means for classifying the relative importance of individual species to 
the structure and stability of a community when experimental manipulations 
are not possible. 

• In chapter 5 I explore the effects of time and space on successional patterns in 
a rocky subtidal community. Using log-linear analysis, I show that transition 
probabilities vary significantly across small spatial scales and over yearly time 
intervals. The implications of these findings for predicting equilibrium species 
abundances and understanding biological processes that drive transient dynam¬ 
ics are discussed.. I also use a set of methods introduced by Cohen et al. (1998) 
to characterize regional variability in the temporal and spatial distribution of 
species. This analysis is conducted to contrast how regional variations in the 
distribution of species compares with local variations in successional patterns 
predicted by log-linear analysis. 

• In chapter 6 I describe the development of a nonlinear Markov chain of a rocky 
subtidal community. The model is parameterized using maximum likelihood 
methods to estimate density-dependent transition probabilities from the Gulf 
of Maine data set. I analyze the best fitting Markov chains to study the effects 
of species interactions on community dynamics, to characterize the sensitivity 
of community dynamics to initial conditions, and to identify multiple stable 
states in the subtidal system. I also examine how changes in the strength of 
interactions among species affect the behavior of the model. 

• Finally, Chapter 7 summarizes the main results as they relate to the general 
theme of this thesis, and proposes some further research directions. 
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Chapter 2 


Habitat fragmentation and 
extinction thresholds on fractal 
landscapes 


The earth belongs to the living, not to the dead. 
-Thomas Jefferson 


2.1 Introduction 

A species living in a fragmented landscape, only part of which is suitable for oc¬ 
cupancy, faces two challenges for persistence. It must balance mortality and repro¬ 
duction within a patch to maintain itself locally, and must locate those parts of the 
landscape that are suitable. As landscapes become more fragmented due to human 
disturbance, the second challenge becomes a critical conservation and management 
issue. 

Habitat fragmentation is important because there exists a threshold level of suit¬ 
able habitat, below which the population goes extinct, even though its vital rates 
are capable of supporting positive population growth in the suitable areas. This was 
first shown by Lande (1987) in an analysis of the Northern Spotted Owl, which is 
1 This chapter was published in Ecology Letters (1999) 2:121-127. 
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endangered because logging has fragmented its old-growth habitat. This result has 
become a paradigm in conservation biology (e.g. McKelvey et a 1. 1993; Lamberson 
et a1. 1992). 

Lande’s model is an ingenious application of demographic theory, but does not 
directly describe the dynamics of occupied and unoccupied territories. He incorpo¬ 
rated the fraction of suitable territories as a factor in a survivorship function, and 
calculated the net reproductive rate by substituting this survivorship function into 
Lotka’s integral equation of population growth. By setting the net reproductive rate 
to 1, he solved for the minimum fraction of suitable habitat permitting population 
growth. 

An alternative approach is to use metapopulation (Levins 1969) or patch-occupancy 
(Caswell and Cohen 1991a, b, 1995) models to investigate the affects of habitat frag¬ 
mentation. Such models describe the landscape as a mosaic of discrete patches and 
focus on the balance between colonization and local extinction. Several single-species 
metapopulation models have shown that a population cannot persist when the fraction 
of suitable patches or territories in the landscape falls below a critical threshold (Nee 
and May 1992; Lamerson and Carroll 1993; Kareiva and Wennergren 1995; Noon and 
McKelvey 1996). Nee and May (1992) extended this approach to a two species com¬ 
petition model and found that habitat destruction decreased the frequency of patches 
occupied by the superior competitor, but surprisingly increased the frequency of the 
inferior competitor. Multi-species competition models (Tillman et a 1. 1994; Stone 
1995) have shown that the most vulnerable species to habitat loss are the top com¬ 
petitors (given that they have lower colonization rates), and that species extinction 
may occur decades after suitable territory has been destroyed (the so called extinction 
debt; Tillman et a1. 1994). 

Because these models assume that every patch interacts equally with every other 
patch, the explicit arrangement of patches has no effect on the results, and such 
models tell us nothing about how the spatial arrangement of habitat destruction 
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effects a population. These effects can be important (Doak et a1. 1992; Mckelvey et 
a1. 1993; Bascompte and Sole 1996). 

In this paper, we derive a discrete time patch-occupancy model (PO model) which 
captures the essential features of Lande’s model and the other single-species models in 
the literature. We then embed this model in an explicit spatial arrangement of habitat 
destruction by transforming it into a stochastic cellular automata model (Caswell and 
Etter 1992, 1999). We compare the results of the two models to study the effects of 
the spatial arrangement of habitat destruction. 

2.1.1 Lande’s Model 

Suppose that the landscape is divided into patches the size of an individual territory, 
and that a proportion h of these are suitable and a proportion u are unsuitable. Let 
p denote the proportion of the suitable territories that are occupied. Suppose that 
individuals with a suitable territory have a survivorship function l'(x) and maternity 
function b(x). At some age before maturity, juveniles inhabit their parents territory 
with probability £ or disperse in search of suitable territory. A juvenile can search as 
many as m potential territories before dying. The probability of finding an unoccupied 
and suitable territory is 

1 - (1 — £){u + ph) m (2.1) 

Since dispersal happens before maturity, the survivorship function for ages beyond 
maturity is simply l'{x) multiplied by (2.1). Thus the net reproductive rate is 

TOO 

Rq = 1 — (1 — e){u +ph) m / l'{x)b{x)dx (2.2) 

Ja 

where a is the age of maturity. Denoting the integral term by Rq (i.e. the net 
reproduction rate conditional on finding a territory), the condition for equilibrium is 

1 - (1 — s)(u +ph) m R! 0 — 1 (2.3) 
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The only variable in this equation is p; the value which satisfies it is the equilibrium 

_ kzk 

h (2.4) 

0 

/i-sV 7 ” 

* - (15) 

Since p* — k when h = 1, Lande identifies k as the demographic potential of the 
population. 

Plotting p* as a function of h for different values of k reveals a clear extinction 
threshold; if h is too small the population goes extinct. 

2.1.2 Metapopulation Models 

Levins (1969) envisioned metapopulation dynamics as a tradeoff between the colo¬ 
nization rate C of empty patches and the local extinction rate m of occupied patches. 
If we let x denote the fraction of empty patches, then the instantaneous change in 
the proportion of occupied patches is 

— = Cx( 1 — x) — mx (2.6) 

The non-trivial steady state solution is x* = 1 — *£, so the population persists only if 
C > m. Note that this solution implies that for biologically realistic values of C and 
m the metapopulation is unable to occupied all available patches. 

Habitat destruction can be incorporated into (2.6) by introducing a term D de¬ 
noting the proportion of patches that have been destroyed. In a landscape where not 
all of the habitat is suitable, the change in the frequency of occupied patches becomes 

— = Cx( 1 — D — x) — mx (2.7) 


habitat occupancy 


V = < 


Where 
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where the non-trivial solution is now x* = 1 — D — ^. Thus for the population to 
persist, D must be less than the critical destruction level D c = 1 — Equation 
(2.7) can be solved for the proportion of suitable territories occupied by making the 
substitutions h = 1 — D and p = |, yielding 

^ = Chp( 1 - p) - mp. (2.8) 

At equilibrium p* = 1 — and the population goes extinct when h < 

2.2 The Patch-Occupancy Model 

Lande’s model takes a population-centered view of the landscape, and calculates the 
fraction of occupied territory from consideration of the survival and reproduction 
of individuals. Metapopulation models take a patch-centered view where occupied 
and unoccupied patches change states in continuous time, depending on colonization 
and mortality rates. A patch-occupancy model also takes a patch-centered view, 
however, patches are now described as state variables which change states according 
to a discrete time nonlinear Markov process. 

Consider a landscape composed of patches, each of which can be in one of three 
states 


State 

Description 

1 

Suitable habitat; unoccupied 

2 

Suitable habitat; occupied 

3 

Unsuitable habitat 


Denote the number of patches in state i by n i} and the proportion of patches in 
state i by x,, for i — 1,2,3. In the model, the proportion h of suitable habitat is 
constant; h = xi+X 2 - We denote the proportion of suitable patches that are actually 
occupied by p = x^jh. This is the variable of interest in the conservation context. 
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No one expects to find owls in parking lots; what is important is how much of their 
remaining potential habitat is actually occupied. 

We will speak of the colonization, growth, and local extinction of populations 
within patches, but the model can also be interpreted in terms of the immigration 
and mortality of individuals on territories. 

A population in an occupied territory goes extinct in a unit of time with proba¬ 
bility 6, and survives with probability 1 — <5 (Figure 2.1). Local extinction produces 
an empty patch. The colonization of such a patch is described by a Poisson process. 
Colonists arrive at the empty patch from other occupied patches in the system. Thus 
the probability that an unoccupied territory is colonized is 

C = P[at least 1 colonist arrives] (2.9) 

= \-e~ M 

where M is the mean number of colonist arriving at a territory. That number is given 
as 


,, dn 2 

M = --- 

ni+n 2 + n 3 

= bx 2 

where b is the number of dispersing propagules produced by each of the n 2 occupied 
patches, and rii + n 2 + n 2 is the number of patches over which those propagules are 
distributed. Thus b can be interpreted as a measure of propagule production; we will 
call it the dispersal parameter. In terms of b , the colonization probability can be 
written as 

C = 1 - e~ bX2 (2.11) 


( 2 . 10 ) 
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1-C 1-5 



Unsuitable habitat 


Figure 2.1: Transition diagram for the patch-occupancy model. The circles represent 
the diff erent, possible states of a patch, and the arrows show the possible state transi¬ 
tions. C is the probability that a suitable unoccupied patch is colonized and S is the 
probability that occupied patch goes extinct. See text for details. 
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The dynamics of the system are given by 


Xj(t + 1) = (1 - C)xi (t) +5x 2 (t) 

x 2 (t + l) — Cx-i(t) + (1 - S)x 2 (t) (2.12) 

z 3 (t + 1) = X 2 (t) 

Because the proportion of suitable patches, h, is constant and equals x\ + x 2 , (2.12) 
can be reduced to a single equation describing the dynamics of occupied patches 

x 2 (t + 1) = C(h - x 2 (t)) + (1 - S)x 2 (t) (2.13) 


where C = 1 - e~ bx2{t) . 

Equation (2.13) can be rewritten in terms of the occupancy p of suitable patches 

as 

Pit + 1) = (l - e-W)) (1 - p(t)) + (1 - 5)p(t). (2.14) 

As h decreases the colonization rate decreases and fewer empty patches are colonized 
during each time step. 


2.2.1 Stability Analysis and Extinction Threshold 


Equation (2.14) has two fixed points: a trivial equilibrium at p = 0 and another 
satisfying 


1 - e~ hb P 
1 - e~ hb P + 5 


(2.15) 


The species can invade the system when the trivial equilibrium is unstable. Denote 
the right hand side of equation (2.14) by f[p,h). The boundary between stability 
and instability of the trivial equilibrium is given by 


df(p, h) 


dp 


p =o 


(2.16) 
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The value of h satisfying this equation is h c , the critical proportion of suitable habitat. 
It satisfies 

' (1 + h c b ) -6 = 1 (2.17) 


the solution to which is 

K = l (2.18) 

Populations with a higher dispersal rate or a smaller extinction rate can tolerate more 
habitat destruction. 

A transcritical bifurcation occurs at h= h c , at which the trivial equilibrium (p = 
0) and the nontrivial equilibrium (p) collide and exchange stability. The conditions 
for such a bifurcation are (Wiggins 1990) 


bpe~ hbp (l-p) =0 (2.19) 

e~ hbp [b{ 1 + p(hbp -hb- 2))] ^ 0 (2.20) 

hbe~ hbp [Kb{p - 1) - 2] ^ 0 (2.21) 

where all the derivatives are evaluated at p = 0 and h = h c . It is easy to confirm 
upon substitution that all three conditions are met. 

The nontrivial equilibrium p is always stable when it exists. To see this, note that 
the second derivative (2.21) is negative for 0 < p < 1, so the one-dimensional map 
of equation (2.14) is concave downward (Fig. 2.2). Since the map is always concave 
downward, p can lose stability only by a flip bifurcation, i.e. when df/dp < —1. 
However, this would require that 


d£ 

dh 

d 2 f 

dpdh 

d lL 

dp 2 



5-1 

1 + hb(l -p)' 


( 2 . 22 ) 


But the left side of (2.22) is always positive and the right side is always negative. 
Hence, the inequality never holds, and flip bifurcations are impossible. 
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Suitable habitat (h) 

Figure 2.3: Equilibrium frequency of occupied patches, p, as a function of h for (a.) 
the PO model, and ( b .) the CA (dashed line) when patches are randomly destroyed. 
CA results are averages from five replicate simulations. The error bars represent 95% 
confidence intervals. The curves are for b — 0.2 and 5 = 0.01, 0.05, 0.10, and 0.15. 

Thus, the long-term dynamics of the model are limited to either extinction (when 
h < h c ), or convergence to a stable equilibrium (when h > h c ). The equilibrium 
frequency, p, of occupied patches can be found by solving (2.15) numerically. Figure 
2.3a shows p as a function of h] the pattern agrees with the results of Lande’s model 
for the Northern Spotted Owl (Lande 1988). 


2.3 Adding Landscape Structure 

To study the influence of landscape structure on dynamics we converted the patch 
occupancy model into a cellular automaton (CA) (see Caswell and Etter 1993, 1999, 
Etter and Caswell 1994, for the relationship between patch-occupancy and CA models 
of this sort). In the CA, both the state and the location of patches are followed 
explicitly through time. Each patch evolves in time following (2.12), but interacts 
only with patches in a local neighborhood. 


39 






We simulated the CA on a 256 x 256 grid with periodic boundary conditions. Each 
patch is in one of three states; unoccupied but suitable for colonization, occupied, 
or unsuitable. Patches change state according to a nonlinear Markov chain whose 
transition matrix is 


A = 


1 -C 5 0 

C 1-5 0 
0 0 1 


(2.23) 


The elements of A correspond to the transition probabilities in Figure 2.1. The 
colonization probability C is given by 


C = 1 — exp(—bx 2 ^) (2.24) 

where x^ is the frequency of occupied patches in the neighborhood of the patch under 
consideration. We report results here for a 7 x 7 neighborhood, which is equivalent to a 
dispersal radius of 3 patches in all directions from an occupied patch. The parameters 
in the CA model are the same as those in the PO model: 5 is the extinction probability 
of a local population within a patch, and d is the dispersal parameter. 


2.3.1 Fractal Landscapes 

In a spatially explicit model the pattern of habitat fragmentation may also affect 
population dynamics. The simplest pattern is a random uncorrelated distribution of 
suitable and unsuitable patches. Real landscapes, however, are often fractal, show¬ 
ing patterns of variance or clumping at different spatial scales (Mandelbrot 1982; 
Krummel et. al. 1987; Milne 1988; Sole and Manrubia 1995). 

We created landscapes with fractal patterns of habitat destruction from fractal 
surfaces whose contours are the traces of fractal Brownian motion. If V (t) represents 
the value of a randomly moving trace at time t , then the change, AV = V(t 2 ) — V{t \), 
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over the increment At = — ti, scales as 


AV oc A t H (2.25) 

where the ruggedness of the trace is controlled by the hurst exponent H (Saupe 
1988). In 3-dimensional space the path of the trajectory produces a fractal surface. 
These surfaces are now familiar from computer graphics books and movies and appear 
remarkably similar to real landscape surfaces (Peitgen and Saupe 1988'). 

The intersection of the surface with a horizontal plane ("flooding” the landscape 
to a specific elevation) creates two sets of patches: those above the plane and those 
below it. By defining the patches below the plane as unsuitable and then progressively 
flooding the landscape to higher elevations, we create fragmented habitats in which 
the frequency of suitable patches declines from 1 to 0, and in which the suitable 
patches are self-similar throughout. 

The edges of these patches have a fractal dimension D given by 

D = 2- H (2.26) 

where 1 < D < 2. As D increases, fragment edges become increasingly rough, and 
suitable patches tend to be isolated in small clusters. Their connectivity is minimized 
when D = 2. On the other hand, as D —* 1, fragment edges become smoother, and 
suitable patches are more likely to be located within large contiguous fragments than 
in small isolated clusters. 

We used the midpoint displacement algorithm published in Saupe (1988) to pro¬ 
duce random surfaces with D = 1.9, 1.5, and 1.1. We generated five replicate random 
landscapes for each fractal dimension. Examples are shown in Figure 2.4 with 50% 
suitable habitat. An uncorrelated random landscape, created by allocating suitable 
and unsuitable patches randomly and independently, is also shown for comparison 
(Fig. 2.4d). 
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Figure 2.4: Three Browian fractal landscapes (a-c) produced using the midpoint 
displacement algorithm (Barnsley et al.). Landscapes are depicted with 50% suitable 
habitat. White corresponds to regions consisting of suitable patches. Black areas 
represent destroyed patches. A landscape in which 50% of the patches were randomly 
destroyed (d) is included for comparison. The fractal landscapes were produced using 
the same integer value to seed the random number generator. 
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2.4 Simulation Methods 


Before the start of each simulation, a' proportion of patches in the landscape were 
designated as unsuitable. The location of unsuitable patches were assigned either 
randomly and independently (Fig. 2.4d), or by superimposing the CA onto one of the 
fractal surfaces (Fig. 2.4a-c). The CA was initialized by setting all suitable patches 
to state 2 (occupied), and run until the proportion of occupied patches converged to 
a stable equilibrium. A simulation was considered to have reached equilibrium when 
the absolute value of the slope of a regression line fit to the last 100 values of p(t) 
was less than 0.01. 

A value of b = 0.2 was used for all CA simulations (different values of d give 
qualitatively similar results for the same value of 5/d). Local extinction probabilities 
were varied between S = 0.01 and 8 — 0.19 so that 0 < 8/b < 1. Five replicate 
simulations were performed for each value of h and for each set of parameter values. 

2.5 Results 

Figure 2.3b shows equilibrium frequencies (p) as a function of h for different values 
of 8 on a random uncorrelated landscape. The behavior of the spatial model on such 
a landscape is similar to that of the non-spatial model. As the proportion of suitable 
habitat decreases, p decreases until a threshold is reached at which the population 
goes extinct. The equilibrium frequency p decreases faster in the CA than in the 
PO model, and the critical proportion of suitable habitat required for persistence is 
larger, especially for higher disturbance probabilities (8 > 0.05), but the spatially 
explicit model creates no qualitatively new results. 

When habitat is destroyed in a fractal pattern, the behavior of the CA changes 
dramatically (Fig. 2.5). The equilibrium frequency (p) on a fractal landscape is 
much higher than on a random landscape or in the PO model. For a given amount 
of suitable habitat, p varies inversely with the fractal dimension of the landscape, i.e. 
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Figure 2.5: Equilibrium frequency of occupied patches, p, as a function of h, for the 
CA on landscapes with different fractal dimensions. Equilibrium curves for the PO 
model (dashed line), and the CA with random habitat destruction, are included for 
comparison. CA results are averages from five replicate simulations. The error bars 
represent 95% confidence intervals. Parameter values are <5 = 0.125, and b = 0.2. 
Statistically, values of p for D = 1.9,1.5, and 1.1 are significantly different for h < 
0.95. This pattern is qualitatively similar for all values of 5/b < 0.9. 


the population does much better when suitable patches are more contiguous. This 
pattern is qualitatively the same for all values of S/b < 0.90. 

Figure 2.6 shows the extinction thresholds for the different habitat destruction 
regimes as a function of 5/b. On a fractal landscape, the CA population can persist 
with much less suitable habitat than predicted by the PO model (except when 5/b > 
0.925, in which case the population always goes extinct, even if h = 1). Decreasing the 
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Figure 2.6: Extinction threshold, h c , as a function of d/6 (for 6 fixed at 0.2) on random 
and fractal landscapes. For the CA, h c is defined as the value of h at which p w 0.01. 
Each point represents the average estimated value of h c in five replicate landscapes. 
The dashed fine gives the extinction threshold for the PO model. 

fractal dimension of the landscape further decreases the amount of suitable habitat 
required for the population to persist. This is especially true when local extinction 
rates are higher (0.5 < <5/6 < 0.9). 
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2.6 Discussion 


Spatially implicit models, such as Lande’s demographic model and single species 
patch-occupancy models, have provided conservation biologist with valuable insights 
into the problem of habitat fragmentation. Patch-occupancy models, however, are 
mean field approximations to spatially explicit CA models (Caswell and Etter 1993,' 
1999). Comparisons between the two kinds of models provide insights into how the 
spatial arrangement of patches in a fragmented landscape affect population dynamics. 

The amount of habitat loss that a population can tolerate depends on the spa¬ 
tial arrangement of suitable and unsuitable habitat. On an uncorrelated random 
landscape, the population occupies slightly less territory than predicted by the PO 
model, and is more susceptible to global extinction. The magnitude of this differ¬ 
ence increases with the local extinction rate. This result is similar to work published 
by Bascompte and Sole (1996), who found that a population was more vulnerable 
to the effects of random habitat destruction in a spatially explicit metapopulation 
model than in its spatially implicit counterpart. Such findings would suggest that 
non-spatial models underestimate the effects of habitat fragmentation on population 
persistence. 

In the real world, however, habitat destruction is rarely a completely random 
process. Instead, fragmentation tends to produce suitable territories of varying size 
with irregular boundaries that are often fractal in nature (Krummel et. al. 1987; 
Milne 1988; Scheming 1991; Sole and Manrubia 1995). On such a fractal landscape, 
a population can tolerate greater habitat destruction, and has a lower extinction 
threshold, than predicted by the PO model. Persistence becomes even more favored 
as the fractal dimension of the landscape decreases. As D —* 1, suitable patches 
become more clumped together and the boundaries between suitable and unsuitable 
territory become smoother. This arrangement insures that on a local scale, patch 
dynamics remain relatively unaffected by habitat destruction outside the periphery 
of a cluster of suitable patches, since few propagules produced in these clusters end up 
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in uns uitable territory. Thus a fractal arrangement of suitable patches increases the 
ability of a population to tolerate habitat destruction by facilitating the recolonization 
of empty suitable territory. This conclusion supports previous studies on reserve 
design for territorial species, which show that the likelihood of population persistence 
increases when suitable territories are clustered together (Doak 1989; Carroll and 
Lamberson 1993; Lamberson et al. 1994). 

Our findings are consistent with those of Dythan (1995), who compared different 
habitat loss regimes in a cellular automata counterpart of a spatially implicit compe¬ 
tition model proposed by Nee and May (1992). In their model, Nee and May showed 
that increased habitat destruction progressively decreases the absolute abundance of 
the superior competitor while increasing the relative abundance of the inferior com¬ 
petitor. In CA simulations, Dythan found that when destruction of patches was 
non-random, the relationship between the superior and inferior competitors changed 
in a similar way, but both species were able to persist in less habitat than predicted 
by Nee and May’s model. 

Although we have shown that a population on a fractal landscape can persist 
with a relatively small amount of suitable territory, the results presented here should 
not be taken as justification to continue destroying native habitats. Instead they are 
meant to point out the need to think about spatial structure and spatial scale when 
considering conservation strategies. It is unlikely that all the effects of landscape 
structure have been incorporated into our model. For instance, habitat connectivity 
can have adverse effects on populations by facilitating the spread of contagious dis¬ 
eases and increasing predation pressure. Models of this effect suggest that increasing 
territorial connectivity beyond a critical point can drive a population to extinction 
(Hess 1996). The interaction between such mortality rates and landscape structure 
in spatially explicit models will have to be reserved for future studies. We believe it 
is safe to say, however, that landscape structure plays a vital role in mediating the 
effects of habitat fragmentation on persistence. 
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Chapter 3 


The effects of habitat destruction 
in finite landscapes: A 
chain-binomial metapopulation 
model 


Interestingly, according to modem astronomers, space is finite. This is a 
very comforting thought—particularly for people who can never remember 
where they have left things. 

-Woody Allen 

3.1 Introduction 

Many species living in fragmented habitats exist as an assemblage of local popula¬ 
tions occupying distinct patches that are connected, to varying degrees, by dispersal 
(Levins 1969; Hanski and Gilpin 1997; Hanski 1999). Ecologists call such assemblages 
metapopulations (Hanski and Simberloff 1997). To persist, a metapopulation must 
balance the extinction of local populations with the colonization of empty patches 
(Hanski 1999). 

A major threat to metapopulation persistence is the destruction of habitat due 
1 This chapter been submitted to Oikos for publication. 
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to human activities. Habitat destruction affects the balance between colonization 
and extinction rates by reducing the number of suitable patches, i.e. patches that are 
capable of supporting a local population. As the number of suitable patches decreases, 
stochastic fluctuations in colonization and local extinctions render the metapopulation 
increasingly vulnerable to global extinction (Lande et al. 1998). 

A common approach to investigating the effects of habitat destruction on the 
abundance and persistence of a species has been to use metapopulation models (e.g. 
Levins 1969; Lande 1987; Nee and May 1992; Lawton et al. 1994; Nee 1994; Tillman et 
al. 1994; Stone 1995; Gyllenberg and Hanski 1997; Bascompte and Sole 1996, 1998; 
Hill and Caswell 1999). These models have shown that the proportion of suitable 
patches occupied declines as habitat is destroyed, and that a metapopulation cannot 
persist when the fraction of suitable patches in the landscape falls below a critical 
threshold. 

Most of these models assume that the landscape consists of an infinite number of 
patches, and describe the landscape in terms of the proportion of these patches that 
are occupied. While such models may provide a good approximation in landscapes 
consisting of hundreds of patches, they may fail badly when the number of suitable 
patches is small (Nisbet and Gurney 1982; Hanski 1999). 

To study the effects of habitat destruction in small finite landscapes, we developed 
a stochastic model that describes the landscape in terms of the number, rather than 
the proportion of occupied patches. The model specifies probabilities of transition 
from any number, to any other number of occupied patches. The change in the number 
of occupied patches has a binomial distribution that can be specified in terms of 
colonization and extinction probabilities. Similar models have been historically used 
in epidemiology (Bailey 1957), and have recently been used by Klok and De Roos 
(1998) and Richards et al. (1999). Because dynamics are a sequence of binomial 
probabilities, the models are called chain-binomial. We will refer to our model as a 
chain binomial metapopulation (CBM) model. 
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In the CBM model, the state of the metapopulation at time t is described by the 
number, n, of occupied patches. The state space is the set of all integers 0 < n < S, 
where S is the number of suitable patches in the landscape. The output of the model 
is a probability distribution vector x(t), where Xi(t) is the probability that n = i. The 
CBM model is an absorbing Markov chain, where extinction is the absorbing state. 

In this paper, we analyze the CBM model to study the effects of habitat destruc¬ 
tion in small landscapes. We use methods from stochastic process theory to examine 
how quasi-equilibrium patch occupancy frequencies, mean extinction times, and per¬ 
sistence probabilities vary as a function of the number of suitable patches and the 
ability of propagules to search for suitable patches. We also describe a method for 
parameterizing the CBM model from time series data, and then use this method to 
study how the frequency of parasitoid outbreaks affects extinction times of a host 
butterfly population. 


3.2 Model Description 

The landscape contains N patches, S of which are suitable, where S < N. Only suit¬ 
able patches are capable of supporting a local population. Suitable patches become 
empty when the population occupying them goes extinct, while empty patches be¬ 
come occupied after being colonized and developing their own populations. Habitat 
destruction decreases the number of suitable patches within the landscape. 

We model the dynamics of the metapopulation as a two-part process. The first 
part describes extinction within patches, while the second part describes the coloniza¬ 
tion of unoccupied patches. For each of these processes we derive a transition matrix 
whose elements give the probability that i patches will be occupied after a unit of 
time given that j patches are currently occupied. 
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3.2.1 Extinction matrix 


The number of occupied patches decreases when one or more local populations go 
extinct. A local population goes extinct with probability 5 and survives with proba¬ 
bility 1 — 5. If n patches are occupied before the extinction process, then the number 
of patches occupied after extinction is a binomial random variable with parameters 
(n, 5). The probability, ms(i,j), that i patches are occupied after extinction given 
that j patches were occupied before extinction is 


m s (i,j) = l 


< M 

l ) 


(1 — 5y8i 1 if i < j 


otherwise 


(3.1) 


for i, j < S. Extinction is represented by a (S + 1) x (S + 1), upper triangular matrix 
Mj whose (i,j) element is ms(i,j). For example, when S = 3 and 5 = 0.1, 


Ms 


( 1.000 0.100 0.010 0.001 N 

0 0.900 0.180 0.027 

0 0 0.810 0.243 

^ 0 0 0 0.729 j 


(3.2) 


3.2.2 Colonization matrix 

The number of occupied patches increases when one or more empty patches are colo¬ 
nized. An empty patch is colonized with probability P c and remains uncolonized with 
probability 1 — P c . If n patches are occupied before the colonization process, then 
the number of patches occupied after colonization is a binomial random variable with 
parameters (n, P c ). The probability that n increases from j to i is the probability that 
i — j patches are colonized, where i > j. Thus the probability m c {i,j) that i patches 
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axe occupied after colonization given j patches were occupied before colonization is 


m c (i,j) = { 


s-j N 

i-3 ) 


Pt^l-Pc) 3 -* if i>j 


otherwise 


(3.3) 


for i, j < S. Colonization is thus represented by an (S +1) x (S + 1) lower triangular 
matrix, M c , whose (i,j) element is m c (i,j). 

We assume that colonization is a Poisson process, where the probability that an 
unoccupied patch is colonized is 


Pc 


P[> 1 offspring arrives at patch] 
1 - exp(~/(n)) 


(3-4) 


where f(n) is the mean number of offspring arriving at each patch in the landscape, 
which depends on the number of occupied patches producing propagules. If propag- 
ules disperse randomly, the expected number of propagules arriving at a patch is 

(3.5) 

where 0 is the number of dispersing propagules produced by each occupied patch. 
The parameter 0 measures both reproductive output and colonization ability. Col¬ 
onization depends on N because propagules are dispersed over both suitable and 
unsuitable patches. An example of the colonization matrix for b = 0.5, S — 3, N = 6 
is 

( 1.000 0 0 

0 0.847 0 

0 0.147 0.847 

v 0 0.006 0.154 


0 N 

0 

0 

1.000 j 



M c = 


(3.6) 




Equation (3.5) assumes that the dispersal is random. Later in Section 2.4 we 
generalize f(n) to incorporate propagule searching ability. 

3.2.3 CBM Model 

The transition from t to t+l requires both colonization and extinction. The transition 
matrix is 

A = M S M C . (3.7) 

A is a nonnegative column stochastic matrix of dimension (5 + 1) x (5+1), whose 
elements, a,j, give the probability that i patches are occupied at time t + l given that 
j patches were occupied at time t. Changing the order of the matrix multiplication 
(i.e. M c Ms) changes when you sample the metapopulation, after extinction or after 
colonization, but does not substantially change the long-term behavior of the model. 
An example of the matrix A, using the example matrices (3.2) and (3.6) from above, 
is shown below 

( 1.000 0.086 0.009 0.001 N 

0 0.789 0.157 0.027 , . 

A = (3.8) 

0 0.121 0.723 0.243 

v 0 0.005 0.112 0.729 , 

Note that state 1 (global extinction) is absorbing and that the largest probabilities 
he on the diagonal of A. 

Figure 3.1 shows surface plots of A for larger values of 5. When <5 is small (fig. 
3.1a) transition probabilities are highest near the diagonal of A and drop off sharply 
as you move away from the largest value in each column. As <5 —► 1 (fig 3.1b.), 
however, the largest transition probabilities lie increasingly above the diagonal of A 
and their distribution in each column becomes more spread out. 
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Given j patches occupied at time t 



Figure 3.1: Graphic representation of the transition matrix A for the case S = 40, 
N = 100 and (3 = 0.5. The value of each element a it j is represented by the height of 
the bar in column j row i. a. Transition matrix for 5 = 0.05. b. Transition matrix 
for 5 = 0.4. 
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3.2.4 Search Ability 

The model so far assumes that propagules passively accept the conditions (suitable 
or not) of the first patch they sample. Plant seeds may do so, but many species 
produce propagules that can sample many patches in search for unoccupied suitable 
space before dying. We can modify equation (3.5) for f(n) in this case by assuming 
that a propagule that lands in an occupied or unsuitable patch, leaves that patch and 
samples the landscape again until it either finds an empty suitable patch, or dies. 
The search process is assumed to occur on a fast time scale within the colonization 
process. 

On the first try, the mean number of propagules arriving at each patch is / 3n/N. 
Of these, a fraction 

n+(N-S) 

N 

land in occupied or unsuitable patches. Of these, a fraction (p survive (0 < <p < 1), 
and try again. The mean number of propagules arriving at each patch who are on 
their second attempt is then 




^ n + (JV-S) ) 


(3.10) 


Similarly, the mean number of propagules arriving at each patch who are on their 
third attempt is given by 


j3n 

~N 




n+ (N — S)\ fn+(N — S)\ 

- N - Jfy - N - J 


(3.11) 


Thus the total number of propagules arriving at each patch within the colonization 
interval is 


f(n) 


fin A (5n (n + (N — S)\ 

-N +4 "n{ - N -j 

f3n 

N — 4>{n + N — S) 


+ 0 : 


2 pn 

N 


'n + (N — S)' 
N 


+ 


(3.12) 
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Note that (3.5) is a special case of (3.12) with <f> = 0. 

The average number of patches that a propagule visits depends on the survival 
parameter <j> and the number of suitable and unsuitable patches 

FTNumber of searches) = —-— (3.13) 

1 - 4>q 

where q = 1 — < 1. Thus there are fewer searches when n is small compared 

to S because the propagules have a better chance, on each try, of finding an empty 
suitable patch. As more habitat is destroyed (S decreases), the average number of 
searches increases because the probability of landing in an unsuitable patch increases. 


3.2.5 Model Analysis 

An example of a single stochastic realization of the model specified by A is shown 
in figure 3.2. The number of occupied patches fluctuates, and in this example the 
population avoids global extinction for 1000 iterations. The probability distribution 
for the number of occupied patches satisfies 


x(t + 1) = A x(t). 


(3.14) 


where Xi(t), is the probability that i— 1 patches are occupied at time t. Figure fig:C3f3 
shows an example of x(f) as a function of time, starting from an initial condition in 
which all suitable patches are occupied. Note that as t —» oo, the probability that 
the metapopulation goes extinct —> 1. 

The state n = 0 is absorbing; hence the matrix A can be rewritten in the form 


A = 


1 

e \ 

v° 

T J 


(3.15) 


where e is a 1 x S row vector of extinction probabilities. The matrix T describes 
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Figure 3.2: Stochastic realization of the number of occupied patches in a single land¬ 
scape as a function of time. Parameter values are S = 35, N — 100, j3 = 0.5, 5 = 0.1 
and <f> = 0. 
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Figure 3.3: Time series of the distribution of landscapes with n occupied patches 
starting from an initial condition in which all patches are occupied. Parameter values 
are S — 10, N = 20, (3 — 0.5, 5 = 0.1 and = 0. 

transitions among the transient states (i = 2,3, • ■ •, S + 1). From the matrix T we 
can calculate the quasi-stationary distribution, the expected extinction time, and the 
probability that the metapopulation survives to time t. 

The quasi-equilibrium distribution gives the probability distribution of the num¬ 
ber of occupied patches, given that extinction has not yet happened and will not 
happen for a long time (Seneta 1966). Let g* denote the probability that i patches 
are occupied (i = 1,..., S). The distribution q can be calculated from the right and 
left eigenvectors, w and v, corresponding to the largest eigenvalue of T: 


qi = mvi. 


(3,16) 


The probability vector q, whose elements are g*, is normalized appropriately so that 
Hili ~ 1- The expected number of occupied patches, n, in the quasi-equilibrium 
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distribution is then 


s 

n = '£ i <li- (3.17) 

i= 1 

and the variance in n is 


° 2 (n) = • ( 3 - 18 ) 

Dividing the expected number of occupied patches by the number of suitable patches 
gives 

p = n/S (3.19) 


where p is the expected proportion of suitable patches that are occupied, given that 
the metapopulation has not gone extinct. We will refer to p as the quasi-equilibrium 
frequency. 

The expected time to extinction is determined from the fundamental matrix of T, 
given by 

F = (I — T)- 1 (3.20) 

where I is the identity matrix (Kemeny and Snell 1976). The expected extinction 
time given that j patches are currently occupied, r J? is calculated by summing the 
jth column of F 

Tj = ^2 fij (3.21) 

i 

The expected extinction time for an initial state selected at random from the quasi¬ 
stationary distribution is 

s 

T=J2<li r i- ( 3 - 22 ) 

2=1 

We will refer tor as the mean expected extinction time of the metapopulation. 
Consider a metapopulation initially occupying j patches. We denote by lj{x) the 
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probability that it persists to ”age” x. Let (T 1 )^. be the (i,j) entry of T x . Then 


tW = E(T*)«- (3-23) 

i 

The expected probability that the metapopulation persists to time x for an initial 
state selected at random from the quasi-stationary distribution is 

K x ) = Y,<lM x )- ( 3 - 24 ) 

i =i 

We will refer to l(x) as the mean persistence probability to time x. 

3.3 Results: The Effect of Habitat Destruction 

In this section we illustrate the effect of habitat destruction on the quasi-equilibrium 
frequency p and the expected extinction time f. We also show how changing the 
searching ability of propagules, <f>, effects p and f. Finally, we show how the mean 
persistence probability l(x) varies as a function of time and the number of suitable 
patches. 

Figure 3.4 shows how p depends on the number of suitable patches S. Reductions 
in S cause p to decline. At very low values of S, however, p increases, because p is 
the expected patch occupancy given that the metapopulation has not gone extinct. 
Thus in the extreme case, if there is only a single patch and the population is not 
extinct, that patch must be occupied and p = 1. 

Figure 3.5 illustrates how quasi-equilibrium frequencies are affected by the propag- 
ule searching parameter 0. Increasing 0 increases both p and the amount of habitat 
destruction the metapopulation can tolerate (fig. 3.5a). Increasing 0 also decreases 
the variance, cr 2 (n), of the quasi-equilibrium distribution (fig. 3.5b). Thus patch oc¬ 
cupancy frequencies for a species with good searching abilities are less variable than 
a species with poor searching abilities. 
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Figure 3.4: Quasi-equilibrium frequency, p, as a function of the number of suitable 
patches, S. The landscape size is N = 100. Parameter values are /3 = 0.5 and (j> = 0. 









Figure 3.5: Varying the searching ability of propagules. a. The effect of the search 
parameter <p on quasi-equilibrium frequencies, p. b. The effect of the search parame¬ 
ter, 0, on the variance of the quasi-stationary distribution, a{n) 2 . Parameter values 
are (3 = 0.5, <5 = 0.1 and N = 100. 

Figures 3.6 shows the expected time to extinction, t, as a function of the number 
of suitable patches. When f is plotted on a logarithm scale, the slope is a convex 
increasing line (fig 3.6a), hence f increases greater than exponentially with S. Plotting 
f on a linear scale (fig 3.6b.) illustrates the large effect that small changes in S have 
on f. There is a threshold number of patches S at which f increases abruptly from 
small values to values so large as to be effectively infinite. The consequences of this 
effect are especially important when f is small (note, the point where the lines in 
fig. 6b appear to meet the x axis depends on the scale of the y axis). Figure 3.6 
also shows that increasing the search parameter <p increases the expected time to 
extinction. Given 2 species with similar values of (3 and S but different values of O. 
the species with the better searching ability will persist longer for a given value of S, 
and can tolerate a greater loss of suitable habitat. 

Figure 3.7 illustrates how the mean persistence probability to time x , l(x) depends 
on S. In general, l(x) is always a monotonically decreasing function of time, however, 
the rate of decrease is highly sensitive to changes in S. For example when N = 100, 
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Figure 3.6: The effect of habitat destruction on extinction times, a. Expected ex¬ 
tinction time f as a function of the number of suitable patches, S, and the search 
parameter 0. f is plotted on a logarithmic scale b. Same as fig. 6a except f is now 
plotted on a linear scale. Parameter values are {3 — 0.5, 8 = 0.2, and N = 100. 

/3 = 0.5, 8 = 0.2, and 0 — 0.5, a reduction in S form 60 to 50 reduces the chance of 
surviving at least 500 years from 85% to 14%. Thus in terms of persistence, there is 
a threshold range of habitat destruction over which small changes in S dramatically 
affect the probability that the metapopulation persists for a long period of time. 


3.4 Parameter Estimation 

In this section we describe a method for estimating /3 and 8 from time series data. 
We focus on the case where 0 = 0 (i.e. propagules have no searching ability). The 
same basic principles apply for species in which 0 > 0, however, it may not always 
be possible to get an accurate estimate of f3 when 0 is unknown. One way around 
this problem is to specify a value of 0 first. If information is available on the average 
n umb er of patches visited by propagules, then equation (3.13) can be used to calculate 
0. Otherwise, (3 and 8 can be estimated for a range of likely 0 values. 
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3.4.1 Likelihood function 


We want to estimate parameters for real metapopulations. Suppose data are available 
in the form of a three-way contingency table, where the (*, j, k) entry of the table is 
the number of patches making the transition from state i, to fate j, at time k. The 
table layout for a single time interval (i.e. k = 1) is shown below 



0(t) 

U(t) 

0(t + 1 ) 

A 

C 

U(t + 1 ) 

B 

D 


where 0(t ) and U(t) are the number of occupied and unoccupied patches at time t, 
and the letters A, B, C, and D represent the number of patches that fall into each 
of the transition categories, and N = A + B + C + D. 

We can estimate (3 and 5 by maximum likelihood. The probability that a patch 
undergoes a particular transition is obtained by summing together the probabilities 
of all possible events that could produce the transition. For example, the probability 
that an occupied patch becomes unoccupied during one time interval is 


P[0(t+l)\Vtt)\ = <>(1 - Prj + 5P C 5(1 - P c ) t (»ry 2 «l - P„) + ... (3.25) 

■ 5(1 - Pc) . 

- -JZTgp- ( 3 - 26 ) 


where P c — 1 — exp(—bO(t)/N). The first term on the right side of (3.25) is the 
probability the patch was disturbed (i.e. local population went extinct) and not 
recolonized. The second term is the probability that the patch was disturbed, recol¬ 
onized, disturbed, and not recolonized; and so on. The other transition probabilities 
are calculated similarly and are shown in the table below 



0(t) 

U(t) 

0(t + 1 ) 

mam 


fta&i 

1-SP, 

U(t + 1 ) 

jjaagi 

1 -Pc 

hni 

1-6P, 
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Suppose the transition data is collected over k time intervals. The likelihood 
function for f3 and 8 given the observed data is 


k 

L(f3,8 1Data) = 


t =1 


/ 1-8 \ A{t) 

\1 -SP C ) 


(8(1-P c )\ m 
V 1 - 8P C ) 


f P c (l - 5) \ c{t) ( 1-P C 
\1-8P C ) U -8P C 


D(t) 


(3.27) 

where 7 is an arbitrary constant, 9 denotes the unknown parameters (fJ and <5 in 
this case) and S denotes the observed data. Statistical calculations utilize the log- 
likelihood of (3.27) given by 


In L((3, <5|Data) 


k 

£ A(t) ln(l - 8) + B(t) ln(5(l - P e ) + ... 

t=i 


C(t) ln(P c (l - 8)) + D(t) ln(l - P c ) -N ln(l - 5P C ) (3.28) 


The value of (5 and 8 that maximizes (3.28), is the maximum likelihood estimate of 
the model parameters given the observed data (Edwards 1992). 


3.4.2 An Example 

Lei and H an s ki (1997) and Hanski (1999) recorded five years of patch occupancy data 
(1993 - 1997) on the Glanville fritillary butterfly Meltiaea cinxia in the Aland islands 
(SW Finland). The data is in the form of yearly spatial maps identifying the location 
of occupied and unoccupied patches in a landscape of 63 patches (see fig. 12.9 in 
Hanski [1999]). We reorganized this data into a three-way contingency table, which 
is shown in figure 3.8a. 

Figure 3.8b shows how the number of patches occupied by M. cinxia and a spe¬ 
cialist parasitoid wasp Cotesia melitaearum varied as a function of time in the Aland 
landscape. The wasp increases the mortality of M. cinxia by attacking the larvae 
of its butterfly host. Between 1993-1995 there was a declining trend in the number 
of patches occupied by M. cinxia due to a high incidence of C. melitaearum (Lei 
and Hanski 1997). In 1995, however, C. melitaearum became nearly extinct in the 
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1993-1994 



°(t) 

U(t) 

Total 

isiion 

24 

2 

26 

ERSH 

8 

29 

37 

Total 

32 

31 

63 


1995-1996 



O(t) 

U(t) 

Total 

0(t+1) 

9 

6 

15 

nsn 

5 

43 

48 

Total 

14 

49 

63 


1994-1995 



BBS 

EBI 

Total 

m 

10 

D 

14 

EBB 

16 

33 

49 

Total 

26 

37 

63 


1996-1997 



°(t) 

U(t) 

Total 

0(t+1) 

13 

8 

21 

Hsu 

2 

40 

42 

Total 

15 

48 

63 


b. 



Year 


Figure 3.8: Patch transition data for a network of 63 patches in the Aland islands 
(SW Finland), a. A three-way contingency table showing the observed yearly patch 
transitions for the Glanville fritillary butterfly Meltiaea cinxia from 1993 - 1997. b. 
Number of patches occupied by the butterfly M. cinxia and its specialist parasitoid 
Cotesia melitaearum as a function of time. 
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1993 - 1995 

1995 - 1997 

1993 - 1997 

ft 

0.3488 

0.9198 

0.6080 

8 

0.4548 

0.2826 

0.4083 

f 

4.18 yrs. 

2.69 x 10 7 yrs. 

13.61 yrs. 


Table 3.1: Likelihood estimates of ft and 8 based on time series data for the Glanville 
fritillary butterfly Meltiaea cinxia. Each column shows the estimated values of ft and 
8 for a specific set of observation years. The expected extinction times for the models 
specified by each set of parameter values are shown at the bottom of each column 
(see text for details). 

landscape, allowing M. cinxia to recover between 1995-1997 (Hanski 1999). 

To look at the effects of C. melitaearum on M. cinxia , we estimated ft and 8 using 
data for 1993-1995 only (parasite outbreak), for 1995-1997 only (low parasite abun¬ 
dance), and for the entire data set (19.93-1997). Table 3.1 fists parameter estimates 
for each time period. Prom these estimates we calculated three transition matrices, 
A 0 (outbreak years), A n (non-outbreak years) and A a n (all years); and calculated f 
for each matrix (Table 3.1). The expected extinction time in non-outbreak years is 
6 million times larger than that in outbreak years. The expected extinction time for 
the entire data set, on the other hand, is only slightly larger than for the outbreak 
years. This pattern suggests that outbreaks of C. melitaearum devastate the M. cinxia 
metapopulation and that more than 2 years are required, following an outbreak, for 
the recovery of M. cinxia. 

We studied how the frequency of parasite outbreaks affect extinction times by 
calculating the matrix 

A F = A k n -'A, (3.29) 

where A: — 1 is the number of years between outbreaks and F = 1 /k is the outbreak 
frequency. The expected extinction time specified by A F (f F ) is in units of k years. 
To convert to a yearly time scale, t f is multiplied by k 

f Fy = kf F . (3.30) 
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Parasitoid Outbreak Frequency 


Figure 3.9: The expected extinction time of Meltiaea cinxia as a function of the 
frequency of Cotesia melitaearum outbreaks, a. f is plotted on a logarithmic scale b. 
f is plotted on a linear scale. 

A plot of ?F y as a function of outbreak frequency is shown in figure 3.9. The 
extinction time of M. cinxia decreases greater than exponentially as the frequency 
of C. melitaearum increases, and has a threshold at F ~ 0.25 (fig. 3.9b). If the 
frequency of outbreaks occur once every four years (or less), then the likelihood that 
M. cinxia goes extinct in the near future is relatively certain. If outbreaks are spaced 
five years apart, however, then the expected time extinction time jumps abruptly to 
about 3000 years. Again we see a sharp threshold in f, this time as a function of the 
frequency of a naturally occurring disturbance. Note that the curves in fig. 3.9 are 
probably a best case scenario, since we assumed outbreaks are evenly spaced in time 
and last only one year. 

3.5 Discussion 

Metapopulation models often describe the effects of habitat destruction on persis¬ 
tence by assuming that colonization and extinction take place in an infinite landscape. 
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Such models predict that the metapopulation either goes extinct or persists at a sta¬ 
ble equilibrium, depending on the proportion of suitable habitat. These predictions, 
however, may fail badly when applied to landscapes consisting of only a small number 
of patches (Nisbet and Gurney 1982, Hanski 1999). In a finite landscape a metapop¬ 
ulation can go extinct simply because all subpopulations happen to go extinct at the 
same time. This is a metapopulation analog of demographic stochasticity, and has 
been termed colonization-extinction stochasticity by Hanski (1999). Colonization- 
extinction stochasticity becomes increasingly important as the size of the landscape 
and the number of suitable patches becomes small (Nisbet and Gurney 1982). To 
describe the dynamics of a finite landscape, a model must account for all possible 
fates of patches. 

Quasi-equilibrium frequencies are a measure of occupied patch densities wit hin 
suitable habitat given the metapopulation has not gone extinct. They are predicted to 
decline as the number of suitable patches, S, declines. At low values of S, however, p 
increases as S declines. This phenomenon corresponds to situations where a scientists 
studies a metapopulation on a small spatial scale (only a few or even one suitable 
patch), or chooses regions to study based on the presence of the species. The latter 
might be common, since most biologists are reluctant to begin studying a species in 
a place where it is known to be extinct. 

Several models (e.g. Lande 1987; Lamberson et al. 1994; With and King 1999) 
have shown that increasing the number of times migrants can search for suitable 
patches increases p and lowers the extinction threshold (i.e. the value of h at which 
a metapopulation goes extinct). In the CBM model, increasing the search parameter 
<j> not only increases p but also reduces the amount of variability in the landscape 
probability distribution (fig. 3.5b). A reduction in c 2 (n ) reduces the probability that 
a landscape will have zero occupied patches, thereby increasing the expected time 
to metapopulation extinction. This result is consistent with previous findings that 
species with better search abilities (as measure here by 4>) can tolerate more habitat 
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destruction. 

Mean extinction times are a greater than exponentially increasing function of S. 
This relationship means that even small amounts of habitat destruction dramatically 
reduce f (fig 6b). As the number of suitable patches decreases, there is a threshold 
response to habitat destruction in which f decreases abruptly from effectively infinite 
values (on an ecological time scale) to very small values. This threshold is also evident 
in terms of l(x), where the probability of persistence over a particular time interval x 
can also decline drastically for small changes in S (fig. 3.7). In general, f and l(x) are 
much more sensitive to changes in S than p, and thus provide a better assessment of 
the effects of habitat destruction in finite landscapes. Predictions of metapopulation 
viability based on p do not capture the true risk of global extinction in small finite 
landscapes as there is not the same threshold response to changes in S, as there is 
with f and l{x). Large underestimates in metapopulation persistence times can result 
by concentrating only on occupied patch frequencies. 

The parameterization methods presented here are straight forward, and are ap¬ 
plicable to any species with a metapopulation structure confined to a small network of 
patches. Once the CBM model is parameterized, it can used to predict how p, t, and 
l(x) are effected by a variety of processes, such as habitat loss, parasite outbreaks (as 
in the case of M. cinxia), habitat degradation, etc. As a large majority of the worlds 
threatened and endangered species are confined to small habitats, it is important to 
understand how they will respond to natural and human induced disturbances. The 
CBM model provides a simple but useful framework for studying these processes. 
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Chapter 4 

A Markov chain model of a rocky 
subtidal community: succession 
and species interactions in a 
complex assemblage 


It is a magnificent feeling to recognize the unity of complex phenomena 
which appear to be things quite apart from the direct visible truth. 

-Albert Einstein 


4.1 Introduction 

Maxine hard substrate communities have proven to be ideal systems for studying the 
dynamics of multi-species assemblages. They are highly diverse, patchy communities 
that tend to be stable on large spatial scales but are relatively unstable on small 
spatial scales (Jackson 1977; Connell and Keough 1985; Sousa 1985). Experiments in 
these systems have shown the importance of a variety of biotic and abiotic processes, 
including competition (Connell 1961a,b, 1972; Paine 1974,1976; Jackson 1977; Quinn 

1982, Sebens 1986), predator-prey interactions (Menge and Sutherland 1976; Connell 

1983, Witman 1985), mutualism (Vance 1978; Steneck 1982, Witman 1987), distur- 

lr This chapter been submitted to Ecology for publication. 
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bance (Dayton 1971; Connell and Slatyer 1977; Paine and Levin 1981; Witman 1987; 
Jackson 1991; Hughes 1994), and recruitment (Sutherland 1974; Roughgarden et al. 
1988; Gaines and Bertness 1992; Graham and Sebens 1996). Because these factors 
can interact in complex ways (e.g. direct species interactions can be nonlinear; shifts 
in species abundance can have indirect effects on the abundance of other species with 
which they do not interact), it is difficult to determine which ecological processes are 
most important for controlling community structure (Wootton 1993). 

One method of untangling this complexity is to seek key processes or species whose 
loss would lead to large changes in the structure of the community (Bond 1994). A 
common approach to identifying these factors is species removal experiments (Paine 
1974, 1992; Menge et al. 1994), in which the structure of a community is monitored in 
experimental plots following the removal of a species from the system. This approach 
has been successful in studies of intertidal communities where, for example, it has been 
shown that the exclusion of the starfish Pisaster ochraceus results in the competitive 
elimination of several sessile species by the mussel Mytilus califomianus (Paine 1966, 
1974). Species such as P. ochraceus whose removal produces a dramatic effect are 
termed strong reactors (Macarthur 1972; Paine 1980), or ’’keystone” species (Paine 
1966). Unfortunately, in communities where no clear keystone species exists, it can 
be extremely difficult to quantify species interaction strengths and to characterize 
the relative importance of weak vs. strong interactions (Menge and Sutherland 1987; 
Paine 1992; Goldwasser and Roughgarden 1993; Laska and Wootton 1998). In many 
offshore marine habitats, such as the rocky subtidal zone, species removal experiments 
are impractical or difficult to perform. To characterize the processes governing the 
dynamics of such communities new theoretical approaches are required. 

The combination of variability on small spatial scales and (relative) stability on 
larger spatial scales, typical of marine and hard-substrate systems, suggests the use 
of Markov chains as a description of community dynamics. These models have been 
used to characterize the dynamics of terrestrial forests (Waggoner and Stephens 1970; 
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Horn 1975; Acevedo 1982; Runkle 1981; Masaki et al. 1992), plant co mmuni ties 
(Isagi and Nakagoshi 1990; Aaviksoo 1995), and insect assemblages (Usher 1979). 
Their application in marine systems has been limited to a few studies on oceanic 
fisheries (Salia and Erxini 1987; Grant et. al. 1988; Formacion and Salia 1994) and 
a comprehensive analysis of a coral reef community (Tanner et al. 1994, 1996). 

A Markov chain describes a community as a landscape of patches, each of which 
is in one of a number of possible states. The state of a patch is determined by 
the presence of an individual of a given species (or species group). The model is 
based on a transition matrix whose (i,j) entry gives the probability that a patch in 
state j changes to state i in one time step. Markov chains usually converge to an 
equilibrium probability distribution of patch states, and most authors have focused 
on that distribution as a prediction of eventual community composition. There are, 
however, many other analytical tools available for the study of Markov chains that 
have not been widely used (Caswell and Cohen 1991a,b). In addition, sensitivity 
analyses originally developed for studying matrix population models (Caswell 1989) 
can be modified for use in Markov chains to investigate the effects of perturbations 
of the transition probabilities on model behavior. 

In this paper we develop a simple Markov chain to characterize the dynamics 
of epifaunal invertebrate communities living on subtidal rock walls in the Gulf of 
Maine. Data for the model comes from permanent quadrats monitored over an 8 
year period. We apply an array of analytical methods to the model to gain insights 
into the key processes underlying the observed structure of the community. These 
include a similarity analysis for classifying species into functional groups, several 
stochastic process indices (turnover rates, recurrence times, and first passage times) 
for quantifying successional dynamics at small spatial scales, and a set of sensitivity 
analysis for identifying the important factors and key species influencing diversity and 
community stability. While our focus here is on rocky subtidal co mmuni ties, these 
methods are wholly applicable to any community of sessile organisms, such as plant 
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communities, coral reefs, or rocky intertidal communities. 


4.2 Background: Rocky Subtidal Communities 

Much of our knowledge of marine hard substrate communities comes from research 
conducted in the rocky intertidal zone (e.g. Paine 1966; Menge 1976; Underwood and 
Denley 1984; Sousa 1985; Roughgarden et al. 1988). In comparison, the ecology of 
organisms living in rocky subtidal zone is much less well known (Witman and Dayton 
2000). Because subtidal organisms inhabit regions that axe typically subject to heavy 

y 

fishing pressures and other human disturbances (Witman and Sebens 1992; Steneck 
1997; Watling and Norse 1998) more attention needs to be focused on these regions 
if we hope to protect the integrity of subtidal systems. 

The rocky subtidal zone encompasses the hard substrate habitat stretching from 
the intertidal fringe down to approximately 200 meters in depth (Witman and Dayton 
2000). Communities are typically dominated by either algae or sessile invertebrates, 
which can occupy up to 90% of the available rock substrate at any one time (Sebens 
1985; Witman and Dayton 2000). Diversity in subtidal communities is generally 
high, with many species coexisting on the substrate surface (Witman 1996). Factors 
thought to be important in maintaining subtidal diversity include predation (Ayling 
1981; Duggins 1983; Witman and Cooper 1983), disease (Scheibling and Hennigar 
1997) competition (Osman 1977, Sebens 1986), physical disturbance (Dayton et al 
1970, Witman 1987, Witman and Dayton 2000), spatial heterogeneity (Witman 1985), 
recruitment (Smith and Witman 1999), sedimentation, current flow (Genovese 1996, 
Genovese and Witman 1999) and the richness of the biogeographical species pool 
(J.D.Witman, F. Smith and R.J. Etter, unpublished). 

Space is frequently a limiting factor in rocky subtidal communities (Osman 1977; 
Russ 1982; Sebens 1986). Species typically compete for space by colonizing and hold¬ 
ing on to empty space (Sebens 1986; Keough 1983), by eliminating nearby species 
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through overgrowth competition (Russ 1982; Witman 1996), or by interference (Ve- 
limiron and Griffins 1979; Bruno and Witman 1996). In some subtidal communities 
overgrowth competition has been shown to be hierarchical with clear dominant species 
(Sebens 1986), while in others competition is best described as a network in which 
competitive interactions between species are often reciprocal (Buss and Jackson 1979; 
Russ 1982). 

The slope of the rock surface and depth can also affect subtidal community struc¬ 
ture. At shallow depths, horizontal and gently sloping substrates are generally dom¬ 
inated by macroalgae, while vertical rock walls are dominated by epifaunal inver¬ 
tebrates (Witman and Cooper 1983; Sebens 1985; 1986, Witman and Sebens 1988; 
Bruno and Witman 1996; Witman and Grange 1998). This pattern is probably regu¬ 
lated by multiple factors but is undoubtedly related to higher light levels on horizontal 
substrates, which creates a more favorable environment for the growth and survival 
of macroalgae (Witman and Cooper 1983). With increasing depth, the abundance 
of sessile invertebrates increases and the abundance of macroalgae decreases (Vadas 
and Stenck 1988; Witman and Sebens 1988). Thus differences between horizontal 
and vertical rock wall communities are less distinct at depths greater than 30 meters 
(Witman and Dayton 2000). 

The focus of our study is a vertical rock wall community located at approximately 
30 meters depth on Ammen Rock Pinnacle in the Gulf of Maine (Witman and Sebens 
1988; Leichter and Witman 1997). The data for our model was collected over an 
eight-year period (1986-1994). It consists of a series of photographs chronicling the 
spatial distribution of sessile species on the rock wall substrate through time. Ten 
replicate quadrats, positioned randomly along a 20 meter span of rock wall habitat 
were photographed at least yearly with a Nikonos V camera mounted on a quadrapod 
frame (as in Witman 1985). Color prints were made of the high resolution color 
slides to identify the species of five major taxa of epifaunal invertebrates (sponges, 
sea anemones, ascidians, bryozoans, and polychaetes) and a species of coralline algae 
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Model States 

Species Type 

State ID 

Bare Rock 


BR 

Hymedesmia 1 sp. 

Sponge 

HY 1 

Hymedesmia 2 sp. 

Sponge 

HY 2 

Myxilla fimbriata 

Sponge 

MYX 

Mycale lingua 

Sponge 

MYC 

Metridium senile 

Sea anemone 

MET 

Urticina crassicomis 

Sea anemone 

URT 

Aplidium pallidum 

Ascidian 

APL 

Ascidia callosa 

Ascidian 

ASC 

Parasmittina jeffreysi 

Bryozoan 

PAR 

Idmidronea atlantica 

Bryozoan 

IDM 

Crisia ebumea 

Bryozoan 

CRI 

Filograna implexa 

Polychaete 

FIL 

Spirorbis spirorbis 

Polychaete 

SPI 

Coralline Algae 

Encrusting algae 

COR 


Table 4.1: Subtidal species identified in ten quadrats located at 30 meters depth on 
Ammen Rock pinnacle in the Gulf of Maine. Species are identified in the model using 
the state codes in the right-hand column of the table. 


(Fig. 4.1). A total of 14 species were identified in the quadrats (Table 4.1). 

4.3 Model Structure and Analysis 

We modeled the dynamics of the subtidal community using a Markov chain. The 
model is defined by its transition matrix A, whose elements, a^. give the probability 
that a patch (i.e. a small discrete area on the rock substrate) in state j at time t 
changes to state i at time t+1. The matrix A is nonnegative (all ay > 0) and has 
the property that each column sums to 1 (i.e. A is column-stochastic). 

The state of a patch is defined by the species that occupies it. A patch can also 
be empty (bare rock). Thus, since there are 14 species in our data set, the number 
of states in our model is 15. Patch states are identified by the abbreviations given in 
Table 4.1. 

In the model, the dynamics of the community are described in terms of patch tran- 
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Figure 4.1: Photo quadrat showing most of the 14 species of epifaunal invertebrates 
investigated. The thinly encrusting sponge Hymedesmia species 1 (HY1) dominates 
most of the space in this particular quadrat. The orange mounding sponge Myxilla 
fimbriata (MYX) is also prominent. Species abbreviations are as in Table 4.1. Not 
shown are the sea anemones, Urticina crassicomis and Metridium senile, which occur 
in large aggregations, the polychaete Filograna implexa, and the bryozoan, Idmidronea 
atlantica. 
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sition probabilities. Transition probabilities depend on several ecological processes 
affecting species abundances. 

1. The probability that an empty patch becomes occupied is a function of col¬ 
onization. Colonization of a patch occurs either by successful recruitment of 
larvae onto bare substrate or by growth of individuals into unoccupied space. 

2. The probability that an occupied patch becomes empty is a function of distur¬ 
bance. Disturbance within a patch is a result of predation, physical disturbance, 
disease, or any other process that cause an occupied patch to become empty. 

3. The probability that a patch occupied by species j becomes occupied by species 
i is a function of species replacement. Species replacement can occur either 
directly through competitive overgrowth (Sebens 1985) or indirectly when a 
mortality event is followed by colonization within a single time step (Witman 
1987, 1996). 

4. The probability a that patch occupied by species i remains occupied by species 
i is a function of persistence. Species that are resistant to disturbance and 
competitive replacement have high persistence probabilities. 

The proportion of patches in each state gives a description of the species composi¬ 
tion of the community. Let x(t) be a column vector giving the probability distribution 
of patch states at time t. Then the species composition at time 1 4-1 is given by 

x(f + 1) = Ax(f) (4.1) 

In our model the time interval is one year. This interval was chosen because sub¬ 
annual observations of quadrat photos show minimal variation in species composition. 
Thus multiplication of x(t) by A projects the community vector forward one year in 
time. 
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Figure 4.2: Cartoon showing the method used to collect state transition data from 
quadrat photos. The squares represent photos from a single quadrat taken at yearly 
intervals. A lattice of approximately 600 points was placed on each quadrat photo 
and species transitions under each point were followed through time. Points were 
space approximately 1 cm apart and are considered to represent a 1 cm 2 patch on the 
substrate wall. 

The largest eigenvalue of A equals 1. The corresponding eigenvector, w, gives 
the equilibrium distribution of patch states. If A is primitive, the community will 
asymptotically converge (as t —> oo) to w from any initial condition. The zth element 
of w ( Wi ) gives the proportion of the landscape occupied by species i at equilibrium. 

4.3.1 Parameter estimation 

Data to construct the transition matrix were obtained by superimposing a lattice of 
evenly spaced points over the quadrat color prints (30 x 20 cm) and counting patch 
transitions through time (Fig. 4.2). Approximately 600 points (a single point every 
1 cm) were assayed per quadrat. We chose this scale because it was approximately 
equivalent to the size of the smallest organism in our data set. For simplicity we will 
refer to each point as a patch, where the size of a patch is taken to be 1 cm 2 . Since 
individuals of many of the subtidal species are capable of growing much larger than 
1 cm 2 , a single individual can occupy more than one patch. 
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The transition probabilities were estimated by constructing a two way contingency 
table in which the i, jth entry gives the number of patches that were in state j 
at some time t and in state i the following year. We constructed the contingency 
table by pooling the data among all quadrats and all time intervals. The transition 
probabilities, a^, were calculated as 


n. 


O'ij — 


V 


Tin 


(4.2) 


where 7iy is the number of transitions from state j to state i (the i, jth entry of the 
contingency table), and rij is the total number of transitions starting in state j at 
time t (sum of the jth column of the contingency table). 

By pooling the data, we are averaging over small scale spatial variability and 
small scale temporal variability to produce the best single realization of a Markov 
chain for the subtidal community. The questions we are exploring with our model 
concern the expected behavior of the community, given that the present conditions are 
maintained indefinitely. As we will show, analyses of homogenous Markov chains can 
reveal important information about the processes giving rise to observed patterns of 
species abundances. The effects of temporal and spatial variability on model behavior 
are explored in a separate paper (Hill et al. 2000 in prep). 


4.3.2 Identifying Functional Groups 

The number of states in our model is large. To make the Markov chain more tractable, 
we introduce a method for combining species into groups based on the functional 
similarity of species roles within the community. Typically, species are combined into 
functional groups based on their degree of taxonomic relatedness (e.g., Waggoner 
and Stephens 1970, Saila and Erzini 1987, Tanner et al. 1994). The problem with 
this approach, however, is that taxonomically related species may have functionally 
different effects on community dynamics. 
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To form functional groups we used a compression algorithm for Markov chain 
matrices developed by Spears (1998). The algorithm works by combining a pair of 
states i and j together to create a new state i U j that represents the probability of 
being in either state i or state j. Spears shows that the probability of a transition 
from any state k into state iUj is simply the sum of the individual state probabilities 


Q-iUj.k — &i,k ~b dj,k 


(4.3) 


while the probability of a transition from state i U j to any state k is a weighted 
average of the elements and a^, 




rritOk^ mjdkj 
mi + mj 


(4.4) 


where is the sum of all the elements in the ith row of A. 

The functional similarity among species is determined by measuring the distance 
between the rows and columns of A. The distance between the rows and columns 
associated with states i and j is given by 




nriidj^ mjdi'k 


mi + mj 


l a M a i,j\J 


(4.5) 


(Spears 1998). The value of 6ij can be thought of as a measure of the degree of 
functional dissimilarity between states i and j. Pairs of species with low values of 8ij 
are good candidates for combing into functional groups. 

Combining a pair of states reduces the dimension of the transition matrix by one. 
The algorithm for compressing the transition matrix is as follows: 


1. Calculate the distance <5y, i, j = 1,..., S, between all pairs of states. 

2. Find the pair of states i* and j* for which <5*.,-. = min(<%). 

3. Compute a weighted average of columns i* and j* (using Eq. 4.5)and place the 
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result in columns i* and j*. 


4. Add rows i* and j* and place the result in row i*. 

5. Delete column j* and row j* from the matrix. 

This process can be iterated by applying the algorithm to the new matrix A;.^. 
(where Aj*uj- is the transition matrix for the combined state i* U j*). 

Combing states i* and j* results in the loss of some information from the model. 
We can quantify the loss of information by calculating the percent change ($) in the 
equilibrium distribution as 

4> = — ■ x 100 (4.6) 

£ 

where Wj-u,-. is the dominant eigenvector of Aj-uj., and w' is a vector formed by 
setting element i* of w equal to the sum of elements i* and j * and then deleting 
element j*. 

4.3.3 Rates and patterns of succession 

A Markov chain describes two spatial scales—the regional dynamics of the commu¬ 
nity and the local dynamics of patches. Although the community converges to an 
equilibrium distribution of patches (characterized by w), individual patches change 
state continuously through time. 

To investigate rates of successional change at the spatial scale of a patch, we 
calculated mean turnover rates, Smouchowski recurrence times for each state, and 
the mean first passage times for pairs of states. These indices are a well-known part 
of the literature on Markov chains (e.g., Iosifescu 1980); they provide insights into 
rates of species change within patches and patterns of succession (Caswell and Cohen 
1991a, b). 
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• The mean turnover rate describes the probability that a randomly selected patch 
changes state between t and t + 1, and is given by 


Wi ^ ~ a «) ( 4 - 7 ) 

i= 1 

where w z is the ith element of the dominant eigenvector, and (1 — an) is the 
probability that a patch in state i changes states in the time interval from t to 
t + 1 (Caswell and Cohen 1991b). 

• The Smouchowski recurrence time 0* of state i is the time elapsing between a 
patch leaving state i and then returning to it again. Its mean is given by 


0i = 


1 — Wi 
Wi( 1 - an) 


(4.8) 


(Iosifescu 1980). A particularly informative measure is the Smouchowski re¬ 
currence time for bare rock (0 : ). Its value gives the mean time a patch stays 
occupied once it has been colonized. 

• The mean first passage time is the average time it takes for a patch in state j 
to first reach state i. This measure can provide insights into the relative rates 
of succession. Let denote the mean first passage time from state j to state 
i. The matrix T whose elements are is given by 


T = (V dg ) -1 (I — Z + Z dg E) (4.9) 

where V is a matrix whose columns all equal w, Z is a matrix given by 
Z = [I — (A — V)] -1 (I is the identify matrix), E is a matrix of ones, and the 
subscript dg denotes matrices containing the diagonal elements of V and Z 
(Iosifescu 1980). Note that the mean recurrence time for a state is given by the 
diagonal elements T zi of T. The difference between T iz and 0*, however, is that 


91 



Ta can be heavily influenced by patches remaining in state i from time t to t +1 
(Caswell and Cohen 1991b). Thus, is a better estimate of recurrence times 
than Ta. 


4.3.4 Sensitivity analysis 

We used sensitivity analysis to investigate how species abundance would vary in 
response to changes in the elements of A. The effect of such changes provides a 
measure of how important a given transition probability is to the overall composition 
of the community. To quantify this effect, we derive sensitivity formulas for the 
dominant eigenvector of a Markov chain and for a scalar measure of species diversity. 

Eigenvector sensitivities provide information on how the equilibrium distribution 
of patches w changes in response to changes in A. Tanner et al. (1984) attempted 
such an analysis for a Markov chain of a coral reef community, but their calculations 
are flawed because they failed to account for the fact that the column sums of a 
Markov chain transition matrix must sum to one. 

The dominant right eigenvector (here represented as Wi) of A gives the stationary 
community structure distribution. The sensitivity of each element in w x to changes 
in A can be found using Caswell’s (1989, 2000) formula for scaled eigenvectors 


d w : 


r) W 1 s (9 m. 

mlwill ^2 ll w ill V a mj 

ddij mzjCi d&ij 


(4.10) 


If the dominant eigenvector is already scaled so || wq ||= 1, then 


where 


f) W 1 

a llwill 

<9wi 

- nr. \ 

dciij 

ddij 

w i / . Q 

k dciij 

dw x 


(m) 

V TIT 

ddij 

w j 2^ 

\ X Wm 


(4.11) 


(4.12) 


Wj is the j th element of w t , a, is element i of the left eigenvector v m , and A m is 
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the m th eigenvalue. 

The derivatives da k j/daij in (4.10) are determined so that the change in ay is 
compensated by the other entries in column j of A. While many compensation 
patterns are possible (Caswell 2000), here we used proportional compensation 

d a fcj _ ~ a kj ^ 23) 

ddij 1 dij 


in which the change in is distributed over the other the entries in the jth column 
proportional to their value. 

Eigenvector sensitivities provide insights into how changes in transition probabil¬ 
ities affect community structure; however, they can also be cumbersome to interpret. 
A simpler approach, akin to the concept of eigenvalue sensitivities in demographic 
modeling, is to use equation (4.10) to compute the sensitivity of a summary statistic 
describing community structure. 

A s ummar y statistic commonly used in community ecology is the Shannon-Wiener 
diversity index (H). Using this index, the diversity of the stationary distribution can 
be calculated as 


H = —^WilnWi (4-14) 

i 

The sensitivity of H to changes in is derived by taking the derivative of H with 
respect to a l3 


dH 

O&ij 


-]T(lni£7jfe + l) 

k 


dw k 

ddij 


(4.15) 


where dw k jda i3 is the fcth element of dw/daij. Equation 4.15 provides a way of 
characterizing how changes in the transition probabilities affect a scalar measure of 
community diversity. 
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4.3.5 Effects of species removals 

We characterized the functional importance of each species by calculating the effect 
of its removal on community diversity. Diversity, as defined by the Shannon-Wiener 
index H, is a measure of both species richness and the evenness in species abundance 
(Lloyd and Ghelardi 1964). Because species removal automatically affects richness, 
we want a measure that is sensitive only to the evenness component of diversity. 
Evenness J(S ) for a community with S species is defined by the ratio 


J(S) = 


H 

MS) 


(4.16) 


where 0 < J(S) < 1. 


To simulate the removal of species i, we first set row i and column i of A to zero 
and renormalized all the other columns so they sum to one. Call this matrix, with 
species i removed, A* and its equilibrium distribution w (i) . Then the percent change 
in evenness A./j is 




(4.17) 


where Ji(S— 1) = Hi/ ln(S- 1) and Hi is obtained from wIf the removal of species 
i has no effect on the relative abundance of the rest of the community, AJ, = 0. If 
A Ji < 0, then species % has a positive effect on evenness. If A Ji > 0 then the presence 
of species % has a negative effect on diversity. 

Species removal also affects community resilience (i.e. the rate the community 
converges to equilibrium). We measured resilience using Dobrushin’s coefficient of 
ergodicity 


a (A) = - max I Oij ~ a ik \ 

z J ’ k i=i 

(Dobrushin 1956a,b). The coefficient a(A) satisfies the condition 


(4.18) 


a(A) — sup 

XlX 2 


||Axi - Ax 2 ll ) 

ll x i-x 2 || j 


(4.19) 
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where Xi and X 2 are probability vectors of size 5x1 whose elements sum to 1 
(xi 7 ^ x 2 ). The value of a (A) gives the minimum contraction in the distance between 
two initial vectors after multiplication by A. Thus the lower the value of 5(A), the 
faster the community is predicted to converge to equilibrium. 

Dobrushin’s coefficient 5(A) gives the minimum contraction rate, but does not 
show how that minimum relates to the typical contraction rate. To quantify this 
typical rate, we calculated the distribution for the contraction rate 


|jAxi - Ax 2 [ 

t|xi — x 2 | 


(4.20) 


by randomly selecting 10000 pairs of probability vectors, Xi and x 2 , from a uniform 
distribution on the 5 — 1 simplex. We represent the mean of this distribution as 4>. 

If the removal of species i increases 5(A) or 4> then species i has a stabilizing 
affect on the community. If the removal of species i decreases 5(A) or 4> then species 
i has a destabilizing affect on the community. 

The measures 5 (A) and 0 give measures of community residence by estimating 
contraction ratios starting from any initial condition. We chose these measures be¬ 
cause disturbances in ecological systems can be quite large (Dayton 1971; Witman 
1987). 


4.4 Results 

The transition matrix A is shown in Table 4.2. The matrix is organized so that 
colonization probabilities are given in the first column (Bare Rock —> Species i) and 
disturbance probabilities are given in the first row (Species i —» Bare Rock). Per¬ 
sistence probabilities are located along the diagonal of the matrix (bold elements), 
while species replacement probabilities are given by the off diagonal elements of A 
(excluding the first row and column). 

The first column of A shows that colonization probabilities are highly variable, 
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ranging over almost two orders of magnitude. The best colonizer by far is the bry- 
ozoan Crisia (a^ = 0.32), which is three times as likely as the next best colonizer 
(.Hymedesmia 1; ai ;2 = 0.10) to occupy a previously empty patch. The remaining 
species in the co mmuni ty have much lower colonization rates (probabilities range be¬ 
tween 0.005 and 0.062). The probability that an empty patch will not be colonized 
by any species over a one-year period is 0.29. 

Persistence probabilities and disturbance probabilities also vary a great deal among 
species, reflecting a wide range of variation in invertebrate life histories. The species 
with the largest persistence probabilities are Urticina (0.863) and Mycale (0.839), 
while the species with the smallest persistence probability is Spirorbis (0.03). The 
species with the largest disturbance probability is Metridium (0.266), while the species 
with the smallest disturbance probability is Hymedesmia 1 (0.029). 

Approximately 94% of the replacement probabilities (off-diagonal elements of A) 
and 90% of the symmetric pairs of replacement probabilities are non-zero (a i3 and 
a 3i > 0; i,j 7 ^ 1), indicating that competition between species is reciprocal (i.e. 
species i can replace species j, and species j can replace species i). Most replacement 
transitions, however, are improbable, with 84% of them having probabilities less than 
0.05. Replacement probabilities greater than 0.05 are generally much larger than 
their symmetrical pair (a i3 » aji or visa versa), indicating that species competition 
is directional (i.e. species i replaces species j more often than j replaces i ). The one 
exception is for the species Hymedesmia 1 and Crisia, which tend to replace each 
other with approximately equal probability. 

The row sum (RS = J2& i a ij) of the elements of the ith row of A gives a measure 
of the competitive ranking of species i relative to other species in the community. 
The highest ranked competitors are Crisia (RS = 2.21), Hymedesmia 1 (RS = 1.71) 
and Filograna (RS = 0.50). The lowest ranked competitors are Parasmittina (RS = 
0.061), Spirorbis (RS = 0.062), and Urticina (RS = 0.067). 
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Figure 4.3: Predicted and observed species abundances, a. Predicted equilibrium 
abundances of subtidal species (given as the percent cover of substrate) at Ammen 
rock in the Gulf of Maine, as given by the dominate eigenvector w of the transition 
matrix A. b. Observed species abundances at three rock wall sites located nearby the 
quadrats for the years 1987, 1990, and 1992. The observed values are means taken 
over all 3 sites and across all 3 years. Error bars represent 1 S.D. 

4.4.1 Equilibrium distribution 

The equilibrium distribution of patch states is shown in Figure 4.3a. At eq uili bri um 
the model predicts that sponges {Hymedesmia spp., Myxilla , and Mycale ) will occupy 
over 50% of the rock substrate, and that Hymedesmia 1 and the bryozoan Crisia will 
dominate the community (together occupying 58% of the substrate). As a comparison, 
Figure 4.3b shows the mean relative abundance distribution at six quadrats located 
within 75 meters of the model quadrats. Mean abundances are for the years 1987, 
1990, and 1991, and include only those species represented by our model (about 
20% of the substrate was occupied by species not found in the model quadrats). The 
equilibrium abundances predicted by the Markov chain fall within the range of species 
variability observed among the six independent quadrats. 
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4.4.2 Quantifying functional relatedness and combining species 
into functional groups 

To q uant ify the degree of relatedness between species in the subtidal community, we 
generated a dissimilarity dendogram based on Spears’ distance metric 5% (Fig. 4.4). 
The dendogram was constructed using the UPGMA method (unweighted pair group 
method with arithmetic mean) (Ayala 1976). The most similar states are bare rock 
and the polychaete Spirorbis (S = 0.022). This implies that the other species do 
not perceive a patch occupied by Spirorbis differently than an unoccupied patch. The 
next most similar states are the ascidians, Aplidium and Ascidia (<5 = 0.040), followed 
by the pairing of the ascidians with the bryozoan Idmidronea (5 = 0.081). The 
dendogram shows there is much more functional redundancy among the bryozoans, 
ascidians, and polychaetes, then among the sponges and sea anemones, and that 
species do not necessarily cluster according to taxonomic groups. 

Figure 4.5 shows the percent change, <&, in the equilibrium distribution (Eq. 4.6) 
as a function of the reduction in the number of states in the model. We reduced 
the size of the state space by combining pairs of states in order of their similarity 
(<5y), using the compression algorithm described in section (4.3.2). As more states 
are combined $ increases, rising to a maximum of about 20% after combing 13 out 
of the original 15 states (at which point A is a 2 x 2 matrix). 

We set a threshold level of $ = 1% (dashed line in Fig. 4.5) as a cutoff criteria 
for combining species into functional groups. Reducing the number of states by five 
produces a change in the equilibrium distribution of about 0.7%. The order in which 
these states are combined by the compression algorithm is as follows: 

• BR with SPI 

• APL with ASC 

• [APL U ASC] with IDM 
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Figure 4.4: Distance dendogram of species states generated using Spears’ 6 {j metric. 
States are clustered according to their degree of functional similarity with other states 
measured as the distance between the rows and columns of the transition matrix 
A (see text for details). The dendogram was created using the UPGMA method 
(unweighted pair group method with arithmetic mean). 
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Figure 4.5: The percent change in equilibrium distribution ($ tI7 ) as a function of the 
reduction in the number of states in the Markov chain. The number of states were 
reduced by combining states in the order of their similarity, (see text for details). 
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Functional Group (FG) 

APL 

ASC 

IDM 

PAR 

_SPI_ 

Table 4.3: Subtidal species that were combined into a single functional group. These 
species were combined because they perform functional similar roles in co mmu nity 
dynamics (see text for details). Species are identified by the state codes used in table 
Tab:C4Tl. 

• [APL U ASC U IDM] with [BR U SPI] 

• [BR U SPI U APL U ASC U IDM] with PAR 

where [i U j U k] means that the states i, j, and k form a functional group. For 
our model, the compression algorithm reduces the size of the state space by combing 
states into a single functional group. While the threshold of 1% is somewhat arbitrary, 
we felt it was a good cut off point because it kept a polychaete state in our model 
( Filograna ). 

We chose to combine the states SPI, APL, ASC, IDM and PAR into a single 
functional group, thereby reducing the number of states in our model to 11. We 
chose not to include BR in this grouping since it represents processes related to 
disturbance, and is a non-species state. We will identify the functional group species 
by the symbol FG (Table 4.3). The compressed transition matrix obtained by forming 
this functional group is shown in Table 4.4. All results throughout the remainder of 
the paper are based on this matrix (which we will still refer to as A). 

4.4.3 Characterizing patch dynamics 

Turnover rates 

The mean turnover rate of a patch is 0.3762yr -1 (Eq. 4.7). Thus on average, a patch 
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Figure 4.6: Smoluchowski recurrence times, a. Bar graph of the Smoluchowski 
recurrence time 9 { for each state in the model. The height of each bar is the expected 
time (in years) elapsing between a patch leaving states i and then returning to it 
again, b. Smoluchowski recurrence times as a function of the equilibrium distribution 
of patch states w. A least square regression of the data on a log-log scale yields the 
power function 9 X = 2.78w^ im ( R 2 = 0.81 p < 0.001). 

changes state once every 2.6 years and approximately 38% of the patches in a given 
quadrat change state each year. 


Smoluchowski recurrence time 

Figure 4.6a shows that Smoluchowski recurrence times 9{ are highly variable among 
species, ranging from 7.2 to 285.9 years. Hymedesmia 1 and Crisia have the shortest 
recurrence times (~ 8 years). The sea anemone Urticina has by far the longest 
recurrence time, while Metridium, Mycale, Hymedesmia 2, and coralline algae have 
recurrence times on the order of 100 years. Figure 4.6b shows that recurrence times 
are highly correlated with equilibrium abundances. Thus the dominant species have 
the shortest return times once they are displaced from a patch. 

The Smoluchowski recurrence time for bare rock, 9 X , gives an estimation of the 
mean time that a patch remains occupied once it is colonized. The value of 9 X for 
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the subtidal system is approximately 16.5 years (Fig. 4.6a). Since the mean turnover 
time is 2.6 years, this result suggests that after colonization a patch changes state 
approximately 6-7 times before it becomes unoccupied again. 


First Passage Times 

Figure 4.7a depicts the elements of the first passage time matrix T. In general, there 
is little variation in first passage times across the rows of T (except for the diagonal 
elements To), indicating that the expected time it takes for a patch to become occupied 
by state i is independent of by the current state of the patch (given that it is not 
currently in state i ). Figure 4.7b illustrates this point by comparing the average 
value of the i th row of T (excluding t«) to 0*. The values of the two measures are 
almost identical suggesting that successional processes are only slightly affected by 
the current state of the community. 

4.4.4 Sensitivity Analysis 

The eigenvector sensitivities characterize how changes in the elements of the transition 
matrix affect the predicted equilibrium distribution in the community. Applying the 
eigenvector sensitivity formula (Eq. 4.10) to the matrix A gives the sensitivity of 
w to chang es in each of the elements of A. Since there axe 11 states in our model 
the n umb er of eigenvector sensitivities is 11 x 11 x 11 = 1331. Table 4.5 shows the 
distribution in the magnitude of the eigenvector sensitivities. The majority of them 
are small. Only about 4.5% of the sensitivities have magnitudes that are less than 
-0.25 or greater than 0.25, suggesting that changes in the majority of transition 
probabilities have little effect on community composition. 

The sensitivity of diversity better illustrates how changes in the transition proba¬ 
bilities affect community structure. Figure 4.8 shows the sensitivity of H to changes 
in the elements of A. We have separated positive and negative sensitivity values 
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a. b. 


u) 



Figure 4.7: First passage times, a. A 3-D bar graph showing the value of the elements 
of the first passage time matrix T. The height of each bar represents the expected 
time (in years) it takes a patch to go from state j to states i. b. Comparison of the 
mean of the i th row of T (excluding ry) with Error bars represent 95% confidence 
levels of row means. 

onto two graphs for visual clarity. Changes in a %3 that have the largest positive af¬ 
fect on H primarily involve transitions from Hymedesmia 1 and Crisia (Fig. 4.8a). 
Diversity is also positively sensitive to several transitions from Bare Rock (i.e. col¬ 
onization probabilities), the FG functional group and Myxilla. Changes in a i3 that 
have the largest negative affect on H primarily involve transitions to Hymedesmia 1 
and Crisia (Fig. 4.8b). Note that the magnitudes of the negative sensitivities are 
much less than the positive sensitivities. In general, the diversity sensitivities predict 
that increases in the replacement probabilities of Hymedesmia 1 and Crisia by the 
less dominant species (especially Urticina and Mycale ) would result in the greatest 
increase in community diversity. 
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a. 



Figure 4.8: Sensitivity of equilibrium diversity (H) to changes in a^-. (a.) Positive 
sensitivities. Bars with large positive values represent transition probabilities that 
lead to greater community diversity when increased, b. Negative sensitivities. Bars 
with large negative values represent transition probabilities that lead to a reduction 
in community diversity when increased. 
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Sensitivity Range 

Number 

Percent 

< -1 

4 

0.30 

-1 to -0.25 

21 

1.58 

-0.25 to -0.1 

60 

4.51 

-0.1 to 0.1 

1173 

88.13 

0.1 to 0.25 

39 

2.93 

0.25 to 1 

29 

2.18 

> 1 

5 

0.38 

Total 

1331 



Table 4.5: Distribution of eigenvector sensitivities. The first col umn gives the range 
of sensitivity values, the second column is the number of sensitivities that fall within 
a given range, and the third column is the percent of sensitivities that fall within a 
given range. There are a total of 1331 eigenvector sensitivities values (11 x 11 x 11) 
for our model. 


4.4.5 Species removal effects 

Figure 4.9 shows the effect of removing a species on evenness. The removal of Myxilla 
produces the largest decrease in A J 2 , however, the overall percent change in even¬ 
ness is only about 5%. The removal of Crisia or the functional group species (FG) 
also reduces evenness, resulting in a decrease in AJ { of about 3%. The removal of 
Hymedesmia 1 , on the other hand, produces the largest increase in A J i: but again 
the overall change in evenness is only about 5%. Thus the removal of any one species 
appears to have little impact on the evenness of the remaining species in the commu¬ 
nity. 

Figure 4.10 shows the effect of removing a species (or bare rock) on the minimum 
contraction rate d(A) and the mean contraction rate 4>. The dashed lines are the 
values of <5 (A) and when all species are present. The greatest increase in ci (A) 
occurs when bare rock is removed, while the greatest decrease occurs when Urticina 
is removed. The greatest increase in occurs when bare rock and Crisia are removed, 
while the greatest decrease occurs when Urticina and Mycale are removed. Thus the 
model predicts that Urticina and Mycale have a destabilizing effect on the community, 
while bare rock, and Crisia have a stabilizing effect on the co mmuni ty 
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State Removed 


Figure 4.9: Predicted change in species evenness diversity, AJ*, resulting from the 
removal of state i from the community. 
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State removed 


Figure 4.10: Predicted change in community resilience resulting from the removal 
of state i from the community, a. Dobrushin’s coefficient of ergodicity, a (A), as 
a function of the removal of species i. a(A) gives the minimum contraction rate 
between two probability vectors after multiplication by A. b. Mean contraction rate 
for random pairs of probability vectors ^ as a function of the removal of species i. 
Error bars represent 95% confidence intervals. The dashed lines in the graphs are the 
a (A) and <p values for the entire subtidal community. 
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4.5 Discussion 


Markov rhains are powerful tools for analyzing the dynamics of sessile communities. 
This paper has presented an array of analytical methods that ecologists can apply to 
Markov transition matrices to quantify the functional relatedness of species, to orga¬ 
nize species into functional groups, and to characterize rates of successional change. 
We also introduced a method for calculating the sensitivity of the dominant eigen¬ 
vector and species diversity measures to changes in species transition 'probabilities. 
Finally, we presented methods for examining the consequences of species removal on 
co mmuni ty diversity and resilience. 

The transition matrix of the Markov chain provides insights into the complexity of 
interactions within the community. In particular, the replacement probabilities (off- 
diagonal elements of A) provide a measure of species connectivity, sensitive to both 
the number and strength of species interactions. Tanner et al. (1994) have proposed 
that the proportion of positive off-diagonal elements of A is an index of the level of 
connectivity in a community, and have suggested that it could be used to characterize 
and compare different ecosystems. While we agree with their assessment, we suggest 
that a better measurement of connectivity would take into account both the number 
and relative magnitude of the positive off-diagonal elements. One possibility is to use 
a mat rix entropy index, normalized to the number of species states in the model. Let 
A be the transition matrix in which the column and row corresponding to bare rock 
have been removed (we are only interested in species replacement probabilities), with 
the col umns renormalized so they sum to one. The normalized entropy index is 

E = ~ (?§ &ij NHN) (4 ' 21) 

where N is the number of states in A, and jVln(JV)) is the maximum entropy of the 
mat rix. A value of E = 1 means that all the replacement probabilities equal 1/N. 
A value of E = 0 means that in each column of A either the persistence probability, 
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or one of the replacement probabilities equals 1 and all the other elements are zero. 
The value of the entropy index for the uncompressed subtidal transition matrix is 
E = 0.40, and for the compressed matrix is E = 0.36. Entropy values for published 
transition matrices of several marine and terrestrial systems are shown in Table 4.6. 
Calculations of Tanner et al.'s complexity ratio ( RC ) are also shown as a comparison. 
There is a slight correlation between E and RC (r p = 0.52, p « 0.11), however, 
this correlation is not significant. More importantly, entropy values appear to do 
a better job of distinguishing between community types than RC. For example, 
three transition matrices published by Rego et al. (1993), describing the dynamics of 
deciduous forests subject to different fire regimes, all have RC values equal to one, 
yet their E values range from 0.43 to 0.75 (Table 4.6). 

The simplest form of information obtained from a Markov chain is the steady 
state distribution of patches. The dominant eigenvector (w) predicts the expected 
equilibrium abundance of species. In our analysis we found that predicted eq uili brium 
abundances were in good agreement with observed species abundances within a 75 
meter radius of the permanent quadrats. Several ecologists have utilized Markov 
chains primarily as a tool for predicting the future composition of a community (e.g. 
Formacion and Saila 1994; Aaviksoo 1995; Srinath 1996, Childress et al. 1998). In 
some cases these predictions have done a relatively poor job (e.g. Childress et al. 
1998). We feel, however, the real strength of Markov chains for community ecology 
lies not in their ability to accurately predict the future, but in the insights they can 
provide into community processes. 

Identifying functional groups: Spears (1998) distance metric <% is a novel way 
of categorizing species into groups based on their level of functional similarity, where 
function is measured by the species role in community dynamics. When this metric is 
applied to the subtidal transition matrix, the results suggest there is a large amount 
of functional similarity between the polychaete Spirorbis, the ascidians Ascdia and 
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Community Entropy (E) Tanner’s complexity ra tio Matrix size Author _ 

Deciduous Forest 0.38 0.85 5 x 5* Waggoner and Stephens (1970) 

Deciduous Forest 0.63 0.80 11 x 11* Horn (1975) 
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transition matrices that do not include a bare rock (empty space) state. See text for details. 



Aplidium, and the bryozoans Parasmittina and Idmidronea. These species were all 
functionally similar to bare rock suggesting that they are relatively weak competitors 
for space. 

The dendogram in Fig. 4.5 shows that it is not always wise to combine species 
into functional groups based solely on taxonomic considerations. For example, the 
sea anemone Metridium is more functionally similar to the sponges Hymedesmia 2 
and Myxilla then it is to the sea anemone Urticina. This suggests that Metridium and 
Urticina play functional distinct roles within the subtidal communities. Incorporating 
these species into a single functional group would have a large effect on model behavior 
and could potentially obliterate important interactions. 

We were able to reduce the number of states in our model by poo ling similnr 
species into functional groups using Spears’ (1998) compression algorithm. This al¬ 
gorithm allowed us to combine states together with only a minimal loss of ecological 
information (less than a 1% change in the equilibrium community). In the subtidal 
community, the compression algorithm reduced the size of the state space by combin¬ 
ing species into a single functional group. The reduction process, however, need not 
follow this pattern. If similar pairs of species cluster into distinct groupings, then the 
compression algorithm will organize species into multiple functional groups. The size 
and number of functional groups formed depends on the particularities of the species 
in a given community and the threshold amount of error (4>) one is willing to accept 
in the model. 

Rates and patterns of community change: Turnover rates and Smoluchowski 
recurrence times Qi describe rates of succession at the spatial scale of a patch, and 
are useful ways of comparing the time scale of local dynamics between different types 
of communities. For example, we calculated turnover rates for transition matrices 
compiled by Tanner et. al. (1994) of three coral reef communities and obtained values 
of 0.10, 0.20 and 0.23 per year. For our model, we calculated a turnover rate of 0.38 
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suggesting that species turnover in the subtidal community is 1.7 to 3.8 times greater 
than in the reef communities studied by Tanner et al (1994). Similarly, comparisons 
of the rate of recurrence of species within patches could be made by calculating the 
mean value of 0* among all species (9). For instance, in the rocky subtidal community 
6 equals 81.4 ± 23.1 years (1 S.E.), while the coral reef communities have 9 values of 
201.2 ± 64.2 years, 171.4 ± 72.4 years, and 190.8 ± 100.3 years. This suggests that 
the subtidal co mmuni ty of epifaunal invertebrates studied here is more dynamic than 
scleractinian coral communities on the Great Barrier Reef. 

First passage times provide information on patterns of successional change. In 
the subtidal community, the order of species replacement is fairly independent of the 
species currently occupying the patch. In general there is no evidence that coloniza¬ 
tion of patches follows any kind of successional pathway (i.e. species i —* species j —> 
species k ). Thus facilitation does not appear to be an important factor in community 
development. Instead our analysis suggests that patch succession follows an inhibi¬ 
tion model (Connell and Slatyer 1977) in which species colonize a patch and maintain 
it until they die or, in some rare cases, are competitively replaced by another species 
(usually either Hymedesmia 1 or Crisia). 

Sensitivity analysis: The eigenvector sensitivities show that changes in the majority 
of the transition probabilities would have little effect on the structure of the subtidal 
co mmuni ty. While only 59 of the 1331 eigenvector sensitivities had large values, there 
was no obvious pattern relating transition probabilities to changes in community 
structure. The problem with analyzing the eigenvector sensitivities in general is 
that there is often so many of them that their meaning is hard to interpret. The 
eigenvector sensitivity formula for Markov chains, however, is useful for calculating 
s umm ary statistics related to the dominant eigenvector. Here we used this formula 
to derive an equation for the sensitivity of Shannon-Wiener diversity as a function 
of changes in the elements of A. Similarly, we could use equation (4.10) to calculate 
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the sensitivity of several other measures of diversity. For instance, the sensitivity of 
Simpson’s index of diversity (D = 1 - Ei^) is -§£- = -2 

The diversity sensitivities show that changes in transition probabilities associ¬ 
ated with Hymedesmia 1 and Crisia have the greatest effect on community diversity. 
Because these species dominate in the community, increases in the ability of less dom¬ 
inant species to replace them is certain to have a positive impact on diversity. While 
this result is not surprising, other sensitivity predictions axe less obvious. For in¬ 
stance, transitions from Bare Rock are predicted to have a greater effect on diversity 
than transitions to Hymedesmia 1 and Crisia. This pattern suggests that competition 
overgrowth by the two dominant species plays a relatively minor role in structuring 
the community compared to colonization processes. This finding is consistent with 
observations by Witman and Sebens (1990), who found little evidence that vegetative 
encroachment by Hymedesmia spp. had a significant affect on community structure 
in the Gulf of Maine. 

Species removal: The removal of Hymedesmia 1 would produce the largest increase 
in species evenness (AJj = 5%). This is not completely surprising since Hymedesmia 
1 occupies the largest proportion of patches. Figure 4.11a shows the percent change 
in the equilibrium abundance of each state as a result of removing Hymedesmia 1. 
While the abundance of all the remaining species increase, several of the less dominant 
species ( Hymedesmia 2, Mycalle , Urticina, and coralline algae) show greater percent¬ 
age increases than the dominant species Crisia. Still, the overall change in evenness 
is relatively small, suggesting that community composition would not change very 
much. 

The removal of the Myxilla would produce the largest decrease in species evenness, 
but again the impact is small (A Jj = -5%). Figure 4.11b shows that in the absence 
of Myxilla the abundance of the dominant species, Hymedesmia 1 and Crisia , would 
increase by about 10%, while increases in the less dominant species range from 2.5% to 
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a. 


b. 




Figure 4.11: The percent change in species abundance as a as a result of removing a. 
Hymedesmia 1, and b. Myxilla from the model. 

8% (except for Hymedesmia 2) which decreases by about 10%). Overall, the removal 
of any one species produces only a small change in evenness, suggesting that none of 
the species plays a keystone role in maintaining diversity. 

The stability analyses, a( A) and 4>, show that the removal of the Bare Rock 
state and Crisia would result in the largest decrease in the convergence rates. We 
can quantify the relative importance of this decrease on community resilience by 
estimating the time it would take two probability vectors to converge to within a 
distance of 0.01 (i.e. ||xi— X 2 II = 0.01). The maximum convergence time is calculated 
as T max — In 0.01/In a (A), while the mean time is T = In 0.01/ In cj). When all states 
are present, T max = 39.6 years and T — 8.2 years. In the absence of bare rock, 
Tm ax — 56.1 years and T = 10.6 years, which amounts to an increase in T max and T 
of 41.2% and 28.5% respectively. The removal of Crisia results in a 22.1% increase 
in T but only an 8.3% increase in T max . Thus, while bare rock and Crisia are both 
important for community resilience, bare rock appears to be more important in cases 
when perturbations are large. These results suggest that processes that open up space 
(e.g. disease, physical disturbance, predation) and processes that effect colonization 
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(e.g. laxval transport, regional species pool), help to stabilize the community. 

While our focus here is on rocky subtidal communities, these methods are applica¬ 
ble to a wide range of terrestrial and aquatic communities, such as plant communities, 
coral reefs, or rocky intertidal communities. Markov chains axe relatively easy to pa¬ 
rameterize. Data can be collected simply by censusing quadrats at fixed intervals 
and recording the species occupying the space at each of a set of points. Pooling the 
data to estimate transition probabilities averages over small scale spatial and tempo¬ 
ral variability and results in a homogenous Markov chain. While such models ignore 
the effects of temporal-, density-, and spatial-dependence, analyses of homogenous 
Markov chains can provide important insights into successional processes and the 
effects of species interactions on community structure. To assess the importance of 
stochastic and nonlinear effects, non-homogenous Markov models can be developed 
which reflect the temporal-, density-, and spatial-dependent realities of the commu¬ 
nity (Hill et al. in prep.). Whether including these higher-level effects provides new 
insights into the role of species in a community, however, remains to be determined. 
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Chapter 5 


Temporal and spatial variation in 
successional patterns in a rocky 
subtidal community: Statistical 
methods, Markov chain models, 
and Species distributions 


Consistency is contrary to nature, contrary to life. 

-Aldous Huxley 

5.1 Introduction 

Subtidal rock surfaces in New England support highly diverse communities of sessile 
invertebrates that form a tight patchwork of species bound to the hard substrate 
(Sebens 1985; Witman and Dayton 2000). Several processes are thought to be im¬ 
portant in maintaining subtidal diversity, include predation (Ayling 1981; Duggins 
1983; Witman and Cooper 1983), disease (Scheibling and Hennigar 1997) competi¬ 
tion (Osman 1977, Sebens 1986), physical disturbance (Dayton et al 1970, Witman 
1987, Witman and Dayton 2000), recruitment (Smith and Witman 1999), sedimen¬ 
tation (Witman and Dayton 2000), and current flow (Genovese 1996, Genovese and 
1 This chapter been submitted to Ecology Letters for publication. 
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Witman 1999). While the rates of these processes can vary over a range of spatial 
and temporal scales (Witman 1995), there is little information about the effect this 
variation has on successional dynamics and community composition. 

Hill et al. (2000) recently developed a Markov chain model of subtidal community 
dynamics, based on observations of species succession in ten replicate quadrats mon¬ 
itored over an eight-year period. To model the community, we subdivided quadrats 
into an array of small patches, each of which is in one of a number of possible states. 
Species or groups of species were identified with states of the Markov chain, and the 
course of succession within patches was defined by a transition matrix A whose (i,j) 
entry gives the probability that a patch in state j changes to state i in one time step. 

The model developed by Hill et al. (2000) is a homogenous Markov chain. It aver¬ 
ages transition probabilities over small-scale spatial and temporal variability to pro¬ 
duce the best single realization of a Markov chain for the subtidal community. Patch 
transition probabilities, however, characterize ecological processes (e.g., colonization, 
disturbance, competition, predation, persistence) that are affected by changes in en¬ 
vironmental conditions. Because homogenous Markov models ignore the effects of 
spatial and temporal variation on transition probabilities, they are inadequate to 
evaluate the effects of environmental change on successional patterns. 

The transition data collected by Hill et al. can be partitioned into categories 
corresponding to specific locations (i.e. quadrat) and time intervals. If multiple tran¬ 
sition matrices are estimated from counts of observed patch transitions at different 
spatial locations and at different times, then the effects of temporal and spatial vari¬ 
ability on successional processes can be tested using log-linear models (Bishop et al. 
1975, Caswell 1989, 2000). These models are the more powerful descendants of the 
Anderson and Goodman (1957) model, which was original proposed to test the time- 
homogeneity of Markov chains. The Anderson and Goodman (1957) test has been 
applied to community transition matrices (Tanner et al. 1994; Wootton 2000), how¬ 
ever, the log-linear approach is more general and can be extended to more complex 
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experimental designs (Caswell 2000). 

Temporal and spatial variability of the environment not only affects the succes- 
sional dynamics of patches, but also the regional dynamics of the community. We can 
characterize the effects of this variability on regional dynamics by examining how the 
relative abundance of subtidal species varies over time and space. This method was 
proposed by Cohen et al. (1998), and involves constructing a series of probability 
matrices describing the distribution of each species among quadrats through time. 
As far as we know, this method has never been applied to an actual data set. 

In this paper we report on the use of log-linear analysis to test for significant 
effects of temporal and spatial variation on patch transition probabilities in the sub- 
tidal community. We use the information from this analysis as a guide for developing 
sets of Markov transition matrices to predict mean species abundances and to char¬ 
acterize successional dynamics. We also examine the effects of temporal and spatial 
variation on regional dynamics by constructing probability matrices describing the 
distribution of species in space and time. These matrices are analyzed to quantify 
how the spatial distribution of each species varies as a function of time, and how the 
temporal distribution of each species varies as a function of quadrat location (Cohen 
et al. 1998). 

5.2 Methods 

5.2.1 Data collection 

The focus of our study is a community of sessile invertebrates (sponges, sea anemones, 
polychaetes, bryozoans, ascidians) and crustose coralline algae living on rock walls at 
30 - 33 m depth on Ammen Rock Pinnacle in the central Gulf of Maine. The data 
consists of a series of photographs of ten permanent quadrats (each 30 x 20cm 2 in 
area) chronicling the spatial distribution of species and bare rock over a eight year 
period (1986 - 1994). The permanent quadrats were photographed at least annually 
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Model States 

Species Type State Codes 

Bare Rock 

BR 

Hymedesmia 1 sp. 

Sponge HY1 

Hymedesmia 2 sp. 

Sponge HY 2 

Myxilla fimbriata 

Sponge MYX 

Mycale lingua 

Sponge MYC 

Metridium senile 

Sea anemone MET 

Urticina crassicomis 

Sea anemone URT 

Aplidium pallidum 
Ascidia callosa 
Parasmittina jeffreysi 
Idmidronea atlantica 
Spirorbis spirorbis 

Ascidian 

Ascidian 

Bryozoan FG 

Bryozoan 

Polychaete 

Crisia ebumea 

Bryozoan CRI 

Filograna implexa 

Polychaete FI1 

Coralline Algae 

Encrusting algae COR 


Table 5.1: Subtidal species identified in the nine quadrats located at 30 meters depth 
on Ammen Rock pinnacle in the Gulf of Maine. Species are identified in the model 
using the state codes in the left-hand column. 


with a Nikonos V camera equipped with a 15 mm lens and two strobes mounted on a 
rigid camera frame (quadrapod) as described in Witman (1985). Within a sampling 
period, each quadrat was treated as an independent replicate, because individual 
quadrats were separated by a horizontal distance of 1 to 5 meters. 

A total of 14 species, each occupying at least 0.5% of the study area, were recorded 
in the quadrats. We grouped these species into 11 different state categories (table 
5.1). We chose to combine the ascidians, Aplidium and Ascidia , the bryozoans, Paras- 
mittina and Idmidronea, and the polychaete Spirorbis into a single functional group 
based on an objective analysis showing a high level of functional similarity among 
these species. Function similarity is measured by the role these species play in com¬ 
munity dynamics (see Hill et al. 2000 for details). 

Transition data were obtained by superimposing a lattice of evenly spaced points 
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over quadrat photographs and following the state of each point through time. Ap¬ 
proximately 600 points (a single point every 1 cm) were assayed per quadrat. We 
chose this scale because it was approximately equivalent to the size of the smallest 
organism in our data set. For simplicity we will refer to each point as a patch, where 
the size of a patch is taken to be 1 cm 2 . Since individuals from many of the subtidal 
species are capable of growing much larger than 1 cm 2 , a single individual can occupy 
more than one patch. 

5.2.2 Statistical Analysis 

The observed transition data form a four-way contingency table N, in which patches 
are classified by their state ( S ) at time t, their fate (F) at time t+1, the time (T) of 
the observation, and the location (L) of the quadrat. The entry n^ki in cell ( ijkl ) of 
the table gives the number of patches making the transition from state i to fate j at 
time k in location l (Caswell 2000). 

To evaluate the significance of temporal and spatial heterogeneity in the transition 
data, we used log-linear analysis (Bishop et al. 1975, Caswell 1989, 2000). In log- 
linear analysis the logarithm of the cell frequencies of N is modeled as a linear function 
of the effects of F, S, T , L, and their interactions. The significance of the effects of 
time, location, and their interaction are evaluated by comparing the likelihood of 
models including progressively more interactions. 

The null hypothesis in tests for differences among Markov chains is one of condi¬ 
tional independence: given the initial state (S), the fate (F) of a patch is independent 
of time and location (Bishop et al. 1975). Because we consider only hierarchical mod¬ 
els, the null model is denoted FS, STL\ which implies that the cell frequencies are 
modeled as a function of the effects of FS and STL. plus the effects of ST, SL, TL, 
S, T , and L (see the appendix for details). 

The effects of time and location on transition probabilities can be tested in several 
ways as shown in figure 5.1. Beginning with the null hypothesis FS,STL at the top 
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of the figure, we add location effects by including the terms FL and FSL. Their 
effect is evaluated by comparing the models FS,STL and FSL,STL. Similarly, we 
add time effects to the null model by including the terms FT and FST and evaluate 
their effects by comparing the models FS,STL and FST,STL (Caswell 2000). 

We can also test for the effects of location by adding the terms FL and FSL to 
the model that includes the time effects and comparing the models FSL,STL and 
FSL,FST,STL. Similarly, we can test for the effects of time by adding the terms FT 
and FST to the model that includes the location effects and comparing the models 
FST,STL and FSL,FST,STL. Finally, we test for the effects of time and location 
by comparing the models FSL, FST, STL and FSTL (Caswell 2000). 

The goodness-of-fit of a model is measured using the log-likelihood ratio statistic 
G 2 which is asymptotically distributed as x 2 with degrees of freedom equal to the 
difference between the number of cells in N and the number of parameters in the 
model. The goodness-of-fit test compares a model to the saturated model FSLT, 
which fits the data exactly. Tests of a specific interaction are assessed by examining 
the difference in G 2 values (AG 2 ) between models that include and exclude that 
interaction. The likelihood statistic AG 2 is distributed as x 2 with degrees of freedom 
equal to the difference between the degrees of freedom for the two models (Caswell 
2000 ). 

Increasing the number of parameters in a log-linear model will always produce a 
better fit to the data. To determine which model best approximates the mechanisms 
generating the transition data we used the Akaike Information Criteria (AIC)(Akaike 
1973; Burnham and Anderson 1998). For log-linear models AIC is defined as 

AIC = G 2 - 2 (df) (5.1) 

(Christensen 1990) where G 2 is the goodness of fit likelihood statistic and df is the de¬ 
grees of freedom of the model. The log-linear model that minimizes AIC is considered 
the most parsimonious, best-fitting model (Burnham and Anderson 1998). 
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5.3 Results: The effects of time and location on 


patch succession 

A graphic representation of the results of statistical tests for the effects of time and 
location on transition data in the subtidal community is shown in Figure 5.1. The 
effect of location is highly significant whether it is evaluated with the effect of time 
excluded or included in the analysis. Time is also highly significant whether or not 
location is excluded or included in the analysis. The interaction of time and location, 
however, is not significant (P = 0.998). 

Comparison of AIC values shows that the best model is FST, FSL, STL (AIC = 
—7566.2). Thus both location and time have important affects on patch transitions in 
rocky subtidal communities. Incorporation of a time x location interaction, however, 
is unnecessary to explain observed patch transition patterns. Note that the model 
FST, FSL has a much lower AIC value than FST, FST, indicating that location 
has a larger effect on transition probabilities than time. 

5.4 Markov chain models: Species composition 
and community dynamics 

Species composition and successional dynamics in the rocky subtidal com muni ty can 
be studied using Markov chain models. The log-linear analysis suggests that these 
models should incorporate the variability in patch transitions associated with time and 
location, without including a time x location interaction. To explore the implication 
of these results we constructed three sets of transition matrices: 1. a set of 10 spatial 
matrices in which the effects of time are ignored, 2. a set of 8 time-varying matrices in 
which the effects of space are ignored, and 3. a set consisting of a single homogeneous 
matrix in which temporal and spatial effects are ignored. 

Let S represent the set of spatially-varying transition matrices, and let S® denote 
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Figure 5.1: Tests for the effects of time (T) and location (L) in a loglinear analysis of 
the subtidal transition data. Each box designates a particular model, its goodness- 
of-fit G 2 statistic, its degrees of freedom, and its AIC value. The top box gives the 
results for the null model (Fate is dependent only on State). The lower boxes represent 
models that include higher order interactions between Fate and T, L, or both. Terms 
added to each model along with the corresponding changes in G 2 (AG 2 ) and degrees 
of freedom are shown along the arrows. 
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the Ith. matrix in S. The element sfj of S® gives the probability that a patch in state 
j will be in state i the following year, in quadrat l. The transition probabilities for 
are estimated as 

jo Hi n m 

b ij ~ 

Hi Hfe n ijkl 

where is the number of patches that change from state j to state i at time k in 
location l. 

Similarly, let T represent the set of time-varying transition matrices, and let T^) 
denote the Arth matrix in T. The transition probabilities for T (k} are estimated as 




E x, 

l Tlijkl 

Hi Hi n ijkl 


(5.3) 


Finally, let A denote the homogeneous transition matrix. The transition proba¬ 
bilities for A are estimated as 


a ij ~~ 


H k H l n ijkl 

Hf Hi H i nm' 


(5.4) 


5.4.1 Equilibrium Predictions 

To examine what the effects of temporal and spatial variation are on equilibrium 
predictions, we compared: 

• the equilibrium distribution that would result if each of the quadrat-specific 
matrices in the set (<S) were treated as a homogenous Markov chain, with the 
observed species distribution in each quadrat. 

• the equilibrium distribution that would result if each of the time-specific ma¬ 
trices in the set (T) were treated as a homogenous Markov chain, with the 
observed species distribution at each time. 

• the equilibrium distribution for the homogenous Markov chain A with the ob¬ 
served species distribution averaged over time and space. 
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Note that the equilibrium predictions of the different sets of Markov chains depend 
exclusively on the transition probabilities measured in the data set, which are inde¬ 
pendent of observed species abundances in the quadrats (see for example Horn 1975). 
Comparing predicted and observed abundances provides insights into how small-scale 
temporal and spatial variation in transition processes affect the ability of Markov 
rhain models to predict community composition. 

The equilibrium distribution of a Markov chain is calculated using eigenanalysis. 
Let M fc be the kth matrix within the set of spatially-vary (temporally-varying) tran¬ 
sition matrices. If M fc is primitive, its largest eigenvalue equals 1. The corresponding 
eigenvector z k (i.e. the dominant eigenvector) gives the equilibrium distribution of 
patch states, to which the community would converge if the transition probabilities 
specified by remained constant. 

If u/ is the dominant eigenvector of S (i) , then the mean equilibrium distribution 
for the set S is 

q = 4£ u < (5-5) 

L i=i 

where L = 10. Similarly, if v fc is the dominant eigenvector of T (fe \ then the mean 
equilibrium distribution for the set T is 

v = ^ E v * ( 5 - 6 ) 

1 k= 1 

where T = 8. Finally, the equilibrium distribution for homogenous transition matrix 
is given by the dominant eigenvector of A, which we shall denote as w. To measure 
the variance in the equilibrium distributions for the sets <S and T, we calculated the 
variance in each element of u* and v*. 

We can quantify the abilifty of the by comparing The equilibrium predictions of 
the different sets of Markov chains depend exclusively on the transition probabilities 
measured in the data set and not on the initial abundances of species in the quadrats 
(see for example Horn 1975). 
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Figure 5.2 shows comparisons of the equilibrium species distribution predicted by 
the Markov chains vs. relative species abundance within the quadrats. Fig. 5.2a 
compares the mean equilibrium distribution for the set of spatially-varying transition 
matrices (u) with the mean value of the average relative species abundance in each 
quadrat between 1986 and 1994. The agreement is excellent. Figure 5.2b compares 
predicted (uj l = 1,..., L) and observed relative abundances, over all species and all 
quadrats. The correlation is r p = 0.94, p < 0.001. Comparisons of the equilibrium 
distributions for the set of temporally-varying transition matrices (fig. 5.2c, d), and 
the homogenous transition matrix (fig. 5.2e, f) with observed species abundances 
yield similar results. This last result indicates that homogenous Markov chains can 
accurately predict community composition in the rocky subtidal zone. Thus while 
temporal and spatial variation is statistically significant, the effect on predicted abun¬ 
dances is biologically trivial. 

5.4.2 Community Dynamics 

Pooling transition data and eigenvector analyses provide information on long-term 
trends in species abundances, but tell us little about community dynamics. Here, we 
illustrate how temporal variation in transition probabilities affect community dynam¬ 
ics, using the set of time-varying transition matrices to model community dynamics 
as a non-homogenous Markov chain. 

The general form for the non-homogenous Markov chain is given by 

x(t +1) = A t X(t) 

= A t A t _x A t _2 ... A lX (0) (5.7) 

where x(t) is a column vector giving the probability distribution of patch states at 
time t, x(0) is some initial distribution of patch states, and the sequential product 
Ai, A 2 ,.. •, A t represents a single stochastic realization of the environment. 
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Figure 5.2: Predicted vs. observed species frequencies. Quadrat-specific transi¬ 
tion matrices S: a. Equilibrium distribution specified by the set of quadrat-specific 
transition matrices u vs. the mean observed species distribution in each quadrat, 
averaged over space (error bars represent 1 SD). b. Equilibrium distribution specified 
for each quadrat Ufc (k = 1,..., L) as a function of the observed species distribution 
in each quadrat (r p = 0.94, p < 0.001). Time-specific transition matrices T: c. 
Equilibrium distribution specified by the set of time-specific transition matrices v vs. 
the mean observed species distribution in each year, averaged over time. d. Equilib¬ 
rium distribution specified for each time interval v*. (k = 1,..., T) as a function of 
the observed species frequencies at each year (r p = 0.98, p < 0.001). Homogeneous 
transition matrix: e. and /. Equilibrium distribution specified by the dominant 
eigenvalue of A vs. the observed species distribution averaged over time and space 
(r p = 0.99, p < 0.001). 
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The dynamics of the community depend on the sequential ordering of the matrices 
A t . If we assume that the matrices A £ (t = 1,..., oo) are drawn from the set T then 
the sequence of transition matrices in equation 5.7 can be described by a column- 
stochastic transition matrix P, whose (i,j) entry is 

Pa = P (A i+1 = T« | (A t = T<*>) (5.8) 


(Caswell 2000). 

We simulated community dynamics using for two different types of environmental 
variability to specify the sequence of matrices in equation 5.7. 

• Periodic environment; in which the sequence of transition matrices was 
specified by the matrix 


^000 
1 0 0 
0 1 0 
0 0 1 


0 1 N 
0 0 
0 0 
0 0 


\ 0 0 0 ••• 1 0 ) 


(5.9) 


Stochastic environment; in which the sequence of transition matrices were 
chosen independently and with equal probability; i.e. according to the matrix 


P = 


/ i i 
8 8 

1 1 
8 8 


I \ 

8 

1 

8 


i I I ... I # 
\ 8 8 8 / 


(5.10) 


Community dynamics specified by the non-homogeneous Markov chain for the 
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Periodic 


Stochastic 



Figure 5.3: Simulated realizations of the non-homogeneous Markov chain showing the 
transient dynamics of Hymedesmia 1, Crisia , Myxilla , Mycale , and Bare Rock in a 
periodic and stochastic environment. 

five most abundant states, Hymedesmia 1, Crisia , Myxilla , Mycale , and Bare Rock, 
are shown in figure 5.3. Simulations in both environments predict that species abun¬ 
dances are highly variable through time, with frequent switches in dominance between 
Hymedesmia 1 and Crisia (where dominance refers to the species with the greatest 
abundance in the community). 

Figure 5.4 shows frequency distributions of dominance times for Hymedesmia 1 
(fig. 5.4a) and Crisia (fig. 5.4b). Dominance time refers to a continuous length 
of time in which species i ( i = HY 1 or CRI ) is the most abundant species in the 
community. The median dominance time for Hymedesmia 1 is three years, with a 
range of 1 to 19 years. The median dominance time for Crisia is one year, with a 
range of 1 to 7 years. 

Table 5.2 lists temporal correlations between the abundances of several subtidal 
species with Hymedesmia 1 and Crisia. The pattern of correlations indicates that the 
composition of the community changes depending on which species, Hymedesmia 1 
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HY1 CRI 



Dominance Time 


Figure 5.4: Frequency distribution of. the dominance time of Hymedesmia 1 and 
Crisia. Results are based on 10,000 stochastic iterations of the non-homogeneous 
Markov chain using the stochastic environment. 



HY 1 

CRI 

BR 

0.256 

-0.880 

HY 1 


-0.557 

HY2 

0.304 

-0.110* 

MYX 

0.114* 

-0.499 

CRI 

-0.557 


FG 

0.270 

-0.510 

COR 

0.237 

-0.713 


Table 5.2: Product-moment correlation coefficients r p for the temporally-varying 
Markov chain. The second and third columns of the table show r p values for 
Hymedesmia 1 (HYM 1) and Crisia (CRI) vs. the states listed in the first column. 
Correlations are based on 10,000 iterations of the non-homogenous Markov chain in 
the stochastic environment. All correlations axe significant at the p < 0.01 level 
except where marked by a *, in which case p < 0.05. 


or Crisia , is most abundant. Since Hymedesmia 1 and Crisia oscillate out of phase 
with each other (Fig. 5.3), the correlation analysis suggests that species composition 
switches between two points in community phase space. 
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Figure 5.5: Comparison of the long-term temporal average of species abundance for 
the non-homogeneous Markov chain (x) and the equilibrium distribution (w predicted 
by the homogenous Markov chain. Averages for the non-homogeneous Markov chain 
are based on 10,000 iterations in the stochastic environment (the periodic environment 
gives similar results). Error bars represent 1 SD. 

The long-term temporal average of the species composition for the non-homogeneous 
Markov chain can be calculated as 


x = 



(5.11) 


where x*, is a vector whose jth element is the abundance of species j at time k. Figure 
5.5 compares x with the equilibrium distribution, w, for the homogenous Markov 
chain. In general, there is no difference in the predicted mean species abundances 
between these methods of calculation. 

The temporal and spatial variations in transition probabilities, as measured by 
log-linear analysis, reflect the highly dynamic nature of species succession in the 
subtidal co mmuni ty. The effects of this variation on community dynamics cannot be 
understood through eigenanalysis, but instead require stochastic realizations of time- 
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and/or spatial-varying Markov chain to examine the behavior of the community. The 
average species composition of the community (over long time scales), however, can 
be obtained either from the dominant eigenvector of the homogeneous Markov chain 
A or by simulating the non-homogenous Markov chain specified by the matrices T 
and then applying equation 5.11. 


5.5 Temporal and spatial variability in species dis¬ 
tributions 

We have shown that transition probabilities and community dynamics vary signifi¬ 
cantly in time and space, but what about the relative spatial and temporal distri¬ 
bution of individual species? To answer this question we employed a method for 
comparing stochastic matrices introduced by Cohen et al. (1998). The general idea 
is to construct a set of probability matrices Ad, in which the kth matrix of the set 
(M fe ) describes the distribution of species k among locations (time periods) at each 
time period (location). Variations in the spatial (temporal) distribution of species k 
is assesses by measuring the mutual distance between the columns of M*. 

Let Q be a 3 dimensional array, where the ( s,l,t ) element of Q is the number of 
patches containing species s, at location 1, at time t. For each species s, we define the 
LxT column stochastic matrix A s with elements A s (l, t ) = Q(s, l, t)/ Q(s, l, t )), 
and the TxL column stochastic matrix B s with elements B s (t, l ) = Q(s, l, t )/ (J2t Q( s > ^ £))• 
Column t of A s is the probability distribution over locations of species s at time t. 

If the columns of A s are nearly identical then there is little temporal variation in 
the spatial distribution of species s (i.e. relative abundances of species s within the 
quadrats are similar over time). Column l of B s is the probability distribution over 
time of the abundance of species s at location l. If the columns of B s are nearly 
identical then there is little spatial variation in the temporal distribution of species s 
(i.e. relative abundances of species s over time are similar between the quadrats). 
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We can assess the mutual distance between columns in any probability matrix M 
using Dobrushin’s coefficient of ergodicity 

d(M) = ^max | \m : j - m : k\\ p (5-12) 

2 j,k 

(Dobrushin 1956a,b), where m : j is the jth column of M and ||x|| p = (53 |x| p ) 1 ^ p is the 
p norm (p > 1) for the vector x. Calculating d(A s ) (s = 1, ..., 12) gives a measure 
of how much the spatial distribution patterns of species s varies in time. Calculating 
d(B s ) gives a measure of how much the temporal distribution patterns of species s 
varies among locations. 

Figure 5.6 shows a(A s ) and d(B s ) values for all 11 species states. The species with 
the greatest temporal variation in their spatial distribution patterns are Hymedesmia 
2, Metridium, and Filograna (fig. 5.6a). The species with the greatest spatial vari¬ 
ation in their temporal distribution patterns are Hymedesmia 2, Mycale, Metridium , 
Urticina, Crista, Filograna and coralline algae (fig. 5.6b). Note that the spatial vari¬ 
ation in the temporal distribution of Bare Rock is also relatively high. If we assume 
that the frequency of Bare Rock is correlated with levels of predation and on rock 
wall invertebrates, then this result suggest that the temporal frequency of predator 
attacks (and other forms of disturbance) is relatively variable in space. Finally, figure 
5.6 shows that the relative variability in the abundance of Hymedesmia 1, Myxilla, 
and the functional group species (FG) are both temporally and spatially stable. 

5.6 Discussion 

To evaluate the significance of temporal and spatial heterogeneity in patch transition 
processes in a rocky subtidal community, we used log-linear analysis of the observed 
transition frequency tables. These models are the more powerful descendants of 
the original Anderson and Goodman (1957) test, which were introduced to test the 
time-homogeneity of Markov chains. While this test has been applied to matrix 


145 




Species 


Figure 5.6: Temporal and spatial variability of relative species abundances in the 
rocky subtidal community. Results are shown using the 1 norm and the 2 norm 
to calculate Dobrushin’s coefficient of ergodicity (see text for details), a. Temporal 
variation in the spatial distribution of each species. A species with a low value of d(A) 
has a spatial pattern of relative abundance that varies less in time than a species with 
a high value of a (A), b. Spatial variation in the temporal distribution of each species. 
A species with a low value of d(B) has a temporal pattern of relative abundance that 
varies less among quadrats than a species with a high value of a (A). 
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population models (Bierzychudek 1982; Cochran 1986) and community transition 
matrices (Tanner et al. 1994), the log-linear approach is more general and can be 
extended to more complex experimental designs (Caswell 2000). The Anderson and 
Goodman test is a likelihood ratio test for comparing the effects of a single factor (e.g. 
time) on transition probabilities. In essence, using it to investigate a more complex 
experimental design is analogous to using all possible pair wise t-tests in place of 
an ANOVA. While some authors continue to use the Anderson and Goodman test to 
develop co mmu nity Markov chains, it is not an appropriate method for characterizing 
the effects of time and space on transition probabilities, and should be retired. 

A log-linear analysis showed that patch transition probabilities varied significantly 
between locations and between years. Such variation can result from stochastic 
processes such as environmentally induced fluctuations in recruitment, or determin¬ 
istic processes such as density- and spatial-dependent species interaction effects. The 
importance of including temporal and spatial variation in a model of the subtidal 
community, however, depends on the questions being asked. If we are interested in 
predicting long term trends in community composition, over large spatial areas, then 
a homogenous Markov chain (formed by pooling all the transition data across time 
and space) should provide adequate estimates of average species densities. On the 
other hand, if we want to understand the biological processes that drive transient 
dynamics in communities, then we need to develop Markov chains that incorporate 
processes affecting patch transition probabilities over smaller spatial and temporal 
scales. 

In communities with many species it is important to quantify how the distribution 
of each species varies over time and space. The methods of Cohen et al. (1998) 
provide a simple but effective way of characterizing the spatial and temporal stability 
of individual species relative to other organisms in the community. For example in 
the subtidal community, the two dominant species, Hymedesmia 1 and Crisia, show 
somewhat different stability patterns. The distribution of Hymedesmia 1 is highly 
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stable both temporally and spatially, suggesting that fluctuations in its abundance 
are highly synchronized across large spatial scales. The temporal distribution of 
Crisia , on the other hand, tends to vary a lot over small spatial scales, indicating 
that Crisia is more affected by changes in local conditions than Hymedesmia 1. 

Figure 5.6 shows that there is more variability between quadrats in the temporal 
distribution of species, than there is over time in the spatial distribution of species. 
This result suggests location has a greater effect on species distribution than time 
and is consistent with the log-linear analysis, which shows a greater effect of space on 
transition probabilities (excluding FT effects has a lower AIC value than excluding 
FL, Fig. 1). Testing for differences in the variability of species distribution, however, 
is only a qualitative test, and is not as statistically rigorous as log-linear analysis. To 
put Cohen’s approach into practice will require the development of statistical tests. 
Such tests can be developed using analytical methods, resampling techniques and 
randomization procedures (Cohen et. al. 1998). 

One of our goals is to understand the relative extent to which deterministic forces 
are responsible for spatial and temporal fluctuations in the transition probabilities. 
To do this we are currently developing maximum likelihood methods to identify tran¬ 
sitions whose values are dependent on species densities. We are using these methods 
to develop a nonlinear Markov chain model of the subtidal community in order to 
study long term community dynamics, look for multiple stable states and nonlin¬ 
ear dynamics, and examine whether species can be competitively excluded from the 
community (Hill, M., J.D. Witman and H. Caswell, in prep.). 
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Appendix 


In log-linear analysis the log of the cell frequencies in N (n^i) is modeled as a linear 
function of the effects of F, S, T, L, and their interactions. The null hypothesis is 
one of conditional independence: given the initial state, the fate is independent of 
time and location. The resulting model for the null hypothesis is 

log Uijki = U + U F (i) + U S (j) + UT(k) + U L (l) 

+UFS(ij) + UST(ik ) + USL(jl) + UTL(kl) 

+USTL(jkl) (5.13) 

where u is the log of the total number of transitions in the table, us{j) represents the 
effect that the jth initial state has on the cell frequencies in N, u F s(ij) represents the 
effect that the interaction of the jth initial state and the ith fate have on the cell 
frequencies in N, etc. 

To test the effect of location, we add the effect of time to the null model by 
including the terms u FF {ik) and u F sr(ijk) 

log riijki = u + u F (i) + us{j) + UT{k) + ulq.) 

+UFS(ij) + u ST(ik) + U-SL(jl) + U T L{kl) 

+U F T(ik) + U F ST(ijk) + USTL(jkl), (5-14) 

which gives the model FST, STL. To test the effect of time, we add the effect of 
location to the null model by including the terms u F l{u ) and u F sL(iji) 

log rtijki = u + u F (i) + u S (j) + u T {k) + u L (i) 

+U F S(ij) + u ST(ik) + U SL (ji) -i- UTL{kl) 

+U F L(il) + U F SL(ijl) + UsTL(jkl)i (5.15) 
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which gives the model FSL, STL. Finally, to test for the effects of time and location 
we use the model FST , FSL, STL 

log riijki = u + up(i) + us(j) + ux(k) + v>L(i) 

+UFS(ij) + UFT(ik) + UFL(U) + u ST(jk) + u SL(jl) + u TL(kl) 
+UpST(ijk) + UpSUijl) + USTL(jkl) (5.16) 

The parameters for all models are estimated by maximum likelihood methods. See 
Caswell (2000) for a general discussion of this methodology or Silva, Raventos and 
Caswell (1991) for an example in a demographic study. 
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Chapter 6 

A Nonlinear Markov Chain of a 
Rocky Subtidal Community: 
Quantifying the effects of species 
interactions on community 
dynamics 


The community stagnates without the impulse of the individual. The im¬ 
pulse dies away without the sympathy of the community 

-William James 


6.1 Introduction 

Nonlinear species interactions can have important consequences for the structure and 
dynamics of natural communities (Lawton 1992; Paine 1992; McCann et al. 1998). 
Numerous models have shown that nonlinearities in population growth rates and 
interspecific interactions can give rise to chaotic dynamics (e.g. Gilpin 1979; Hastings 
and Powell 1991; Kot et al 1992; Caswell and Neubert 1998), and that community 
stability is dependent on the number and strength of these interactions (May 1973; 
McCann et al. 1998). 

In marine benthic communities, nonlinearities can arise from density-dependent 
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effects on larval settlement (Planes et al. 1998;), post-settlement mortality (Weiss 
1948; Connell 1985; Zajac et al 1989; McShane 1991), and competition for space be¬ 
tween adults (Grant 1977; Denley and Underwood 1979; Sebens 1985,1986; Bertness 
1989; Witman and Dayton 2000). Most empirical studies attempting to quantify the 
effects of species interactions on the dynamics of benthic systems typically involve 
only one or a few sessile species (e.g., Davis 1988, Raimondi 1990, Minchinton and 
Scheibling 1991). Paine (1992) proposed measuring per capita interaction strengths 
among species through a series of species removal experiments. As an alternative, we 
propose to estimate interaction strengths using Markov chain models. 

Markov chains provide a statistical approach to modeling community dynamics 
in marine systems and have been used to characterize successional change in coral 
reef communities (Tanner et al. 1994, 1996), subtidal communities (Hill et al. 2000), 
and intertidal communities (Wootton 2000). A Markov chain describes a community 
as a landscape of patches, each of which is in one of a number of possible states. 
States are defined by the presence of an individual of a given species or species group 
(functional group). The model is based on a transition matrix A, whose (i,j) entry 
gives the probability that a patch in state j changes to state i in one time step. Most 
applications of Markov chains to marine and terrestrial communities have been limited 
to linear models, in which transition probabilities are assumed to remain constant over 
time (e.g. Waggoner and Stephens 1970, Horn 1975, Usher 1979, Acevedo 1982, Salia 
and Erxini 1987, Grant et. al. 1988, Isagi and Nakagoshi 1990, Masaki et al. 1992, 
Rego et al. 1993, Tanner et al. 1994,1996, Formacion and Salia 1994, Hill et al. 
2000 ). 

A Markov chain can incorporate effects of species interactions on patch succession 
by allowing transition probabilities to depend on species densities (Caswell and Cohen 
1991a,b, 1995, Barradas and Cohen 1994, Barradas et al. 1996). If x(t) is a vector 
whose ith element gives the proportion of the landscape occupied by state i at time 
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t, then the dynamics of the community over time is given by 


x(t + 1) = A x x(t) (6.1) 

where the subscript on A x means that the elements of the transition matrix, a^(x t ), 
are functions of the state densities at time t. These functions must take into account 
the requirement that A x is always nonnegative (i.e. Oij(x t ) > 0; V t > 0) and that 
each col umn of A x always sums to 1 (A x is colinnn-stochastic). 

Nonlinear Markov chain models for simple communities have been studied theoret¬ 
ically by Caswell and Cohen (1991a,b, 1993,1995), and Barradas et al. (1996). Those 
models described various models for competition, succession, predator-mediated co¬ 
existence, and local and regional determination of species richness. However, we know 
of no attempt to apply this theory to data on real communities (Although Cornell 
and Karlson (1997) have invoked some of it to explain patterns of species diversity in 
coral reef communities). 

In this paper we develop a nonlinear Markov chain of a rocky subtidal community 
to study the effects of species interactions on community dynamics. We derive maxi¬ 
mum likelihood methods to estimate density-dependent transition probabilities from 
spatial time series data. The data for our model comes from permanent quadrats 
monitored over an 8-year period of epifaunal invertebrate communities living on sub- 
tidal rock walls in the Gulf of Maine. We study the dynamics of the model using 
n um erical simulation, characterizing the behavior of the community starting from 
a large set of initial conditions. We also parameterize nonlinear Markov chains by 
calculating species densities at different spatial scales, to study the dependence of 
interactions strengths on neighborhood size. Finally, we use bifurcation analysis to 
char acterize how changes in the strength of the interactions among species affect the 
temporal dynamics of the community. 
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Model States 

Species Type 

Bare Rock 


Hymedesmia 1 sp. 

Sponge 

Hymedesmia 2 sp. 

Sponge 

Myxilla fimbriata 

Sponge 

Mycale lingua 

Sponge 

Metridium senile 

Sea anemone 

Urticina crassicomis 

Sea anemone 

Aplidium pallidum 

Ascidian 

Ascidia callosa 

Ascidian 

Parasmittina jeffreysi 

Bryozoan 

Idmidronea atlantica 

Bryozoan 

Crisia ebumea 

Bryozoan 

Filograna implexa 

Polychaete 

Spirorbis spirorbis 

Polychaete 

Coralline Algae 

Encrusting algae 


Table 6.1: Invertebrate species identified in the photo quadrat data set. These species 
were originally used by Hill et al. (2000) to develop a linear Markov chain of a rocky 
subtidal community. 


6.2 Data collection 

The focus of our study is a vertical rock wall community located at approximately 30 
meters depth on Ammen Rock Pinnacle in the Gulf of Maine (Witman and Sebens 
1988; Leichter and Witman 1997). The data for our model were collected over an 
eight-year period (1986-1994). They consist of a series of photographs chronic lin g the 
spatial distribution of sessile species on the rock wall substrate through time. Ten 
replicate quadrats, positioned randomly along a 20 meter span of rock wall habitat 
were photographed at least yearly with a Nikonos V mounted on a quadrapod camera 
frame (as in Witman 1985; Cayer et al. 1999). Color prints were made of the 
high resolution color slides to identify the species of five major taxa of epifaunal 
invertebrates (sponges, sea anemones, ascidians, bryozoans, and polychaetes). A 
total of 14 species were recorded in the quadrats (Table 6.1). 

This community has already been described by several authors (e.g., Osman 1977, 
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Model States 

State Codes 

Hymedesmia 1 sp. 

HYM 

Myocilla fimbriata 

MYX 

Mycale lingua 

MYC 

Metridium senile 
Urticina crassicomis 

SEA 

Crisia ebumea 

CRI 

Filograna implexa 

FIL 

Aplidium pallidum 
Ascidia callosa 
Parasmittina jeffreysi 
Idmidronea atlantica 
Spirorbis spirorbis 

FG 

Bare Rock 

BR 


Table 6.2: Species groups used to develop the nonlinear Markov chain. Species groups 
axe identified in the model using the state codes in the right-hand column of the table. 


Ayling 1981, Russ 1982, Duggins 1983, Sebens 1985,1986, Witman 1985,1987,1996, 
Scheibling and Hennigar 1997, Witman and Dayton 2000) and modeled as a linear 
Markov rhain by Hill et al. (2000a). In the linear model five species were combined 
into a single group ( Aplidium , Ascidia , Parasmittina, Idmidronea, and Spirorbis ) as 
they were shown to perform functionally similar roles in community dynamics (see 
Hill et al. 2000 for details). In this paper we reduced the state space further by 
e liminating two species ( Hymedesmia 2 and coralline algae) from the data set (due to 
their low abundance) and combining the sea anemones, Metridium senile and Urticina 
crassicomis, into a single group. Thus the state space for the nonlinear Markov chain 
consists of 7 species states, plus a bare rock state representing empty patches (table 
6 . 2 ). 

Data for the analysis were obtained by superimposing a lattice of evenly spaced 
points over quadrat photographs and following state transitions at each point through 
time. Approximately 600 points (a single point every 1 cm) were assayed per quadrat. 
We chose this scale because it was approximately equivalent to the size of the smallest 
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HYM 

MXY 

From 

MYC SEA 

CRI FIL FG BR 

To: HYM 
MXY 
MYC 

SEA 

On(x t ) 

0 2 i(x t ) 

o 3 i(x t ) 

Oi 2 (x*) 

o 22 (x t ) 
032 (x t ) 

a 13( x t) 

0 2 3 (Xt) 

033 (X t ) 

. Ois(x t ) 

. o 2 s(Xf) 

. o 3 s(x t ) 

CRI 


• 

I 

: 

FIL 

j 

I 

• 

*. : 

FG 

BR 

flsi(xt) 

OS2 (X t ) 

OS3(x t ) ••• 

. o ss (x t ) 


Table 6.3: Schematic representation of the transition matrix for the nonlinear Markov 
chain. The vectors x t in the matrix indicate that the transition probabilities are 
functions of the states densities at time t. 


organism in our data set. For simplicity we will refer to each point as a patch, where 
the size of a patch is taken to be 1 cm 2 . Since individuals of many of the subtidal 
species are capable of growing much larger than 1 cm 2 , a single individual can occupy 
more than one patch. 

6.3 A nonlinear Markov chain model 

A no nlin ear Markov chain is based on a matrix of density-dependent transition prob¬ 
abilities. The matrix is arranged so that column j represents transitions from state 
j at time t to states i, i = 1,..., S, at time t -f 1. The form of the transition matrix 
for the subtidal community is shown in table 6.3. 

Estimating the dependence of transition probabilities on species densities requires 
a flexible functional form and a procedure to estimate the parameters in the function. 
Multiple logistic regression is a familiar method for describing a binary outcome as a 
function of a set of independent variables. Here, we use the generalization of logistic 
regression to polychotomous variables describing multiple outcomes (e.g. Cox 1970, 
Hosmer and Lemeshow 1989). 
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Suppose there axe multiple outcomes for some event with probabilities q t , i = 
1,..., 5, and associated with each outcome is a vector z of independent variables and 
a vector of parameters fy. Polychotomous logistic regression defines a logit function 
for each of the possible outcomes. To ensure that the resulting set of probabilities 
s um to one, we write logit functions by specifying one of the outcomes as a reference 
value. If we let qs be the reference outcome, then the logit function for g* is 


In — = & 0> i + bi,i z i + • • • — +bs t iZs 
Qs 

= (bi,z) 


( 6 . 2 ) 


where the first element (z 0 ) of z = 1 and (b i? z) represents the scalar product of the 
two vectors. The parameter b 0 is a constant, which is unaffected by changes in the 
independent variables. Using the fact that J2iQi = 1, it follows that the probability 
of outcome i is 


Qi = S 


exp((bj,z)) 


1 i ex P(( b i> z )) 


l «?((*>*>*» 


if i - 1,..., S - 1 
if i — S 


(6.3) 


To develop a no nlin ear Markov chain, we treat the probabilities of all transitions 
from state j (i.e column j of A*) as a polychotomous logistic regression problem. 
The vector of independent variables is z = [l,x] r , where x is the vector of species 
densities and the probability qi is the matrix entry g^(x). Since we are interested 
in the effects of species interactions on transition probabilities, we use the bare rock 
state, BR, as the reference outcome. 

Let b^ be the vector of parameters associated with the probability of changing 
form state j th to state i in a single time step. The set of functions a i:? (x) in the jth 
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column of A axe given by 


Oy(x) 


I 


exp 




■>S—1 


i+£ fc _ex p 


1 


i+Efc _ex p 


(( b “' z )) 


if i < S 
if i = S 


(6.4) 


Note that the logistic function can range from nearly linear (no effect of density) 
to having a sharp threshold in the density effects. It is also flexible enough to rep¬ 
resent both positive and negative density effects. Positive effects might result from 
gregarious settlement or from local reproduction; negative effects from competitive 
interactions. 


6.4 Methods 


6.4.1 Parameter estimation 


The transition probabilities are estimated from data on the number of patches that 
change from state j at time t to state i at time t + 1, and the densities of the various 
states in each quadrat at time t. The density of state i within a quadrat is the number 
of patches occupied by state i divided by the total number of quadrat patches. 

Let Bb) be a matrix associated with the jth column of A, whose i th row (denoted 
Bp"*) is the vector b^. The elements of Bb) are estimated by maximum likelihood 
methods. The likelihood function is 


L( Bb)|Data) = mb in 


T L 

nn 

i= 1 1=1 


X 


exp((B?\z tt )) 



/ 

f 1 ^ 

N s ,i,t' 

1 k 

\l + £f 1 exp((B^ ) ,z u))j 

S ' 

J 


(6.5) 


where T is the number of time intervals, L is the number of quadrats, z i t is a vector of 
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species densities in quadrat l at time t, and N iy i yt is the number of patches that change 
from state j to state i in quadrat l within the interval (t, t+ 1). This function insures 
that the probabilities in each column of A x always sum to one, and all possible pairs 
of probabilities have a logistic relationship (Cox 1970). 

The log-likelihood is 

lnL(B^IData) = EE N 1U { B?' } - $ lt ) + ■■■ 
t= 1 z=i 

+N lS -i)u(B<s-i) - *«) - ( 6 . 6 ) 

where $i t = ln{l + Y,k ex P ((B** z it))}- The values of the coefficients in that 
maximize (6.6), are the maximum likelihood estimates of the model parameters in 
col umn j. This process is repeated for each column of A x . 

6.4.2 Set of candidate models 

The set of models we fit to the data represent different hypotheses regarding the 
number of species interactions affecting the replacement of species j by species i. 
These include a model with no species interactions (M 0 ), plus a set of models in 
which each (x) in column j is a function of the density of a single species (Mi), 
of two species (M 2 ), of three species (M 3 ), etc. The number of parameters estimated 
for each col umn of A x depends on which model we are fitting to the data (table 6.4). 
The maximum number of parameters is (S + 1) x (S — 1). This corresponds to the 
situation in which each (x) in column j is a function of the densities of all the 
species in the community (model Mg). 

For a given model, Mh h = 0,... , 8, we used the MATLAB routine fminu to find 
parameter values that maximized the likelihood function (Eq. 6.6). To determine 
which model best approximates the mechanisms generating the transition data we 
used the Akaike Information Criteria {AIC) (Akaike 1973; Burham and Anderson 
1998). If Lh is the likelihood of model M h and pn is the number of parameters in 
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Model Name 

No. of interactions 

No. Parameters 

M 0 

0 

5-1 

Mi 

1 

2(S-1) 

m 2 

2 

2(S-1) 

M s 

S 

(5 +1)(5 — 1) 


Table 6.4: The models fitted to the subtidal data, with the number of species inter¬ 
actions affecting each entry , i = 1,..., 5, in the j-th column and the number of 
parameters estimated for each column of A x . The variable 5 is the n umb er of states 
in the Markov chain (5 = 8). 


model Mh, then the AIC for Mh is 


AIC = — 2 In Lh + 2ph 


(6.7) 


The model that minimizes AIC is considered the most parsimonious, best-fitting 
model (Akaike 1973). Because models are fit separately to each column, it is possible 
to have different best fitting models for different columns of A x . 

For each model Mh, h = 0,..., 8, we assume that the set of the probability 
functions Gy (x) in column j depend on the densities of exactly h species (table 6.4). 
Each ay(x), however, can depend on a different set of species (as long as the set 
contains h species). This means that the total number of possible models for each 
column is 



SI 

n!(5 — n)! 


s 


( 6 . 8 ) 


and for the entire transition matrix is 


(£( 51 
\n=0 V n! (S “ «)! 



s 


(6,9) 


Thus for a community with 8 states (S = 8) the total number of possible nonlinear 
Markov chains is approximately 10 199 . Obviously, this is an incredible large model 
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space to search for the best-fitting nonlinear Markov Chain. Not surprisingly, trying 
to find the best-fitting model using polychotomous logistic regression methods alone 
takes an unreasonably long time. 


6.4.3 Model Fitting using a Two Step Method 

To more efficiently search this model space, we developed nonlinear Markov chains 
using the following two-step process. 

1. We treat each as a binary logistic regression problem to rank the strength 
of the effect each species has on the replacement of species j by species i. 

2. We then use the species rankings obtained in step 1 to specify a set of poly¬ 
chotomous logistic models for each column of A x . 

We can write a binary logistic model for each a^, by assuming there are two 
outcomes for a patch in state j; it either becomes occupied by state i at time t + 1, 
or it does not. The logit function for is 


1*737- = < c ’y> 


( 6 . 10 ) 


where c is a vector of parameters and y = [1, x] is a vector of species densities. Solving 
(6.10) for ctij yields 

p<c,y> 

Oil = 7 7 (6.11) 

3 1 + e< c ’y> 


To rank the species according to their effect on a^, we first assume that a^- depends 
on the density of a single species. For each species k (k = 1,..., S) we maximize the 
likelihood formula 


lnL(c|Data) = EE n (*>^) ((c,yft) - ln(l + e (c ’ yit> ) (6.12) 

i t 
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where y i t axe the densities of species k in quadrat l at time t and n(i,l,t ) is the 
number of patches that change from state j to state i in quadrat l within the time 
interval (t,t + 1). The species that yields the largest likelihood value (In L) is ranked 
as having the strongest effect on the replacement of species j by species i. We denote 
this species as x^ 1 ). 

To identify the seconded ranked species we repeated this procedure, assuming that 
a,ij depends on the density of x^> and a second species m (m ^ x^). The species m 
that yields the largest likelihood value (in combination with species x^ 1 )) is ranked as 
having the second strongest effect on the replacement of species j by species i. This 
process is continued until all species have been ranked according to their effect on 
dij (x). Note that estimations of the parameter values c axe important for determining 
the best ranking of species, however, these values are not used to parameterize the 
nonlinear Markov chain. 

Once the species ranking axe specified for each a^, we write the transition proba¬ 
bilities in each column of A* as a set of polychotomous logistic functions (Eq. 6.4). 
To fit a specific model M h to column j, we create a set of density matrices in the 
form 


z it = 


(1 

1 

1 

H 

4 1} 

... r (b 
x s -1 

U<‘> 

T (W 

x 2 

... r {h) 
X S -1 


(6.13) 


where x* is the density of the species with the Arch strongest effect on the replacement 
of species j by species i (as specified by the binary logistic analysis). For example, to 
fit the model M 2 we substitute the set of matrices 


z a 


T (D T « 

. t ( 2) r (2) 
y x 2 


T (i) 

X S -1 

r (2) , 


(6.14) 
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into the polychotomous likelihood function (Eq. 6.6). 


6.4.4 Best fitting model 

Table 6.5 shows log-likelihood and relative A AIC value (A AIC = AIC value relative 
to the best model) for the models M h , h = 0,..., 8, for each column of A x (HYM, 
MYX, MYC, etc.). The best fitting model for each column has a A AIC value equal to 
0.0. In general, the best fitting models include the densities of all but one or two of the 
species in the co mmuni ty. For example, the best fitting model for the set of functions 
a 1>t -(x) describing the replacement of Hymedesmia by species i (i = 1,..., 8) include 
the densities of seven species. 

Table 6.5 shows that the best nonlinear Markov chain is comprised of a transition 
matrix A x whose first column consists of a set of probability functions specified by 
model My, whose second column consists of a set of probability functions specified 
by model M 6 model, etc. This transition matrix can be represented using the symbol 
M 76777766 , where the first entry in the subscript gives the number of species interac¬ 
tions for column one, the second entry gives the number of species interactions for 
column two, and so on. The rank ordering of species interaction strengths for each 
of the probability functions a i7 (x) is shown in the appendix. 

6.5 Results: Nonlinear Dynamics of the Subtidal 
Community 

6.5.1 Numerical simulations 

The analysis of the nonlinear Markov Chain A /76777766 depends on numerical simula¬ 
tion, as there is no hope of finding analytical solutions for a model of this complexity. 
N um erical analysis were carried out by specifying an initial probability vector x(0) 
(where x(0) is an 8 x 1 non-negative column vector whose elements sum to 1) and 
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Model 

n 

HYM 

logL A AIC 

MYX 

log I AAIC 

MYC 

logL A AIC 

SEA 

logL A AIC 

M 0 

7 

-8739.9 

597.4 

-3213.3 

273.9 

-1457.6 

407.9 

-1681.9 

517.4 

Mi 

14 

-8565.5 

262.6 

-3125.7 

105.3 

-1361.2 

229.1 

-1649.8 

467.1 

m 2 

21 

-8469.4 

84.5 

-3084.2 

43.5 

-1287.2 

95.1 

-1515.5 

4212.6 

M 3 

28 

-8438.5 

36.7 

-3059.9 

9.0 

-1253.1 

41.0 

-1438.4 

72.4 

m 4 

35 

-8428.4 

30.3 

-3051.7 

6.6 

-1237.5 

23.8 

-1413.0 

35.4 

m 5 

42 

-8419.2 

26.0 

-3046.2 

9.6 

-1224.4 

11.6 

-1395.5 

14.5 

Me 

49 

-8411.3 

24.2 

-3034.4 

0.0 

-1217.0 

10.7 

-1385.8 

9.0 

m 7 

56 

-8392.2 

0.0 

-3031.0 

7.2 

-1204.6 

0.0 

-1374.3 

0.0 

m 8 

63 

-8391.3 

12.2 

-3029.1 

10.1 

-1202.8 

10.4 

-1374.2 

13.9 


Model 

n 

CRI 

logL A AIC 

FIL 

logL A AIC 

FG 

logL A AIC 

BR 

logL A AIC 

M 0 

7 

-7862.4 

913.5 

-2456.3 

347.3 

-4435.2 

578.1 

-5432.2 

817.4 

Mi 

14 

-7559.2 

321.0 

-2323.6 

95.9 

-4239.7 

196.1 

-5240.0 

442.9 

m 2 

21 

-7443.8 

104.4 

-2303.9 

70.6 

-4169.4 

74.5 

-5069.7 

120.6 

Ms 

28 

-7418.5 

67.6 

-2269.9 

16.4 

-4145.5 

40.7 

-5038.8 

72.7 

m 4 

35 

-7393.4 

31.6 

-2260.0 

10.6 

-4128.0 

19.6 

-5015.2 

39.5 

M 5 

42 

-7378.5 

15.8 

-2252.1 

8.8 

-4117.9 

13.4 

-5003.3 

29.7 

M 6 

49 

-7371.4 

15.6 

-2244.2 

7.0 

-4104.2 

0.0 

-4981.5 

0.0 

m 7 

56 

-7356.6 

0.0 

-2233.7 

0.0 

-4099.7 

4.9 

-4976.7 

4.5 

m 8 

63 

-7355.5 

11.7 

-2233.6 

13.9 

-4098.6 

16.7 

-4975.4 

15.9 


Table 6.5: The models fitted to the photo quadrat data, with the number of parame¬ 
ters (n) for each column, the log likelihood (log L), and the AIC value relative to the 
best model ( AAIC ). 


then iterating the model according to equation (6.1). At each iteration the value of 
the elements of A x were calculated as a function of x(f) using formula (6.4). 

The nonlinear Markov chain has several classes of potential behavior 

• Convergence to a stable equilibrium distribution x 

• Convergence to multiple stable equilibrium distributions x* * 

• Periodic, quasiperiodic, or chaotic dynamics 

To examine all possible behaviors of the model we ran 10,000 simulations for 500 time 
steps. Each simulation was initialized by randomly choosing a probability vector Xq 
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from a uniform distribution on the 5-1 simplex. 


6.5.2 Equilibrium dynamics 

Simulations of the nonlinear Markov chain revealed three coexisting attractors; a sta¬ 
ble 2-cycle and two stable equilibria (Fig. 6.1). The 2-cycle has the largest basin of 
attraction, with convergence to the 2-cycle attractor occurred in 57.2 % of the simula¬ 
tions (Fig. 6.1a). Species densities on this attractor are almost completely dominated 
by Mycale, which occupies 89.1% of the patches (Fig. 6.1b). The functional group 
(FG) has the second highest density, with an average patch occupancy of 7.1%. Three 
of the species states, Filograna, Crisia and the sea anemones, are almost completely 
eliminated from the community. 

Initial conditions within the second largest basin of attraction converge monoton- 
ically to a stable equilibrium (occurring in 34.2% of the simulations) (Fig. 6.1c,d). 
Species densities on this attractor are dominated by Crisia and Hymedesmia, which 
occupy 55.1% and 30.0% of the substrate, respectively. We refer to this equilibrium 
distribution as the CRI-HYM equilibrium. 

Initial conditions within the smallest basin of attraction converge slowly to a stable 
equilibrium via a damped 2-cycle oscillation (occurring in 8.6% of the simulations) 
(Fig. 6.1e,f). This attractor is dominated by the sea anemones, which occupies 49.6% 
of the substrate at equilibrium (SEA equilibrium). The remaining proportion of the 
substrate is almost completely unoccupied (BR = 49.0%), with no other species 
having densities greater than 0.5%. 

While the 2-cycle has the largest basin of attraction, the species frequencies pre¬ 
dicted by the CRI-HYM equilibrium are much more similar to observed frequencies 
in the subtidal quadrats. Figure 6.2 shows a comparison of the CRI-HYM equilib¬ 
rium species frequencies with the mean observed frequencies. The mean observed 
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Figure 6.1: Coexisting attractors for the best fitting nonlinear Markov chain, in which 
parameters were estimated based on species densities within quadrats. Time series 
for the attractors are shown in a, c, and e, along with the proportion of times the 
model converged to each attractor starting from 10,000 random initial conditions. 
Histograms of the equilibrium species frequencies for each attractor are shown in b, 
d, and /. 
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Figure 6.2: Mean abundance of subtidal species among quadrats over time (Observed) 
vs. CRI-HYM eq uili bri um distribution. Error bars represent 1 SD. The predicted 
CRI - HYM equilibrium abundances of all species except Crisia fall within 1 stan¬ 
dard error of observed abundances. 


frequency of species i was calculated as 


o = Iy-f JttL 


(6.15) 


where n it i ft is the number of patches containing species i in quadrat l at time t. 


6.5.3 Convergence times 

Convergence times for each attractor were calculated by simulating the model M 76777766 
according to Eq. (6.1) starting from 10,000 initial conditions x(0), randomly selected 
from a unif orm distribution on the (5-1) simplex. Convergence times were calculated 
as the mean number of iterations required to satisfy the following criteria 

j2\x i (t-2)-x i (t) |< 0.001 (6.16) 

1=1 
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where Xi(t) is the frequency of state i at time t (we used Xi(t — 2) instead of Xi(t — 1) 
because one of the attractors is a 2-cycle). The mean convergence time for the 2-cycle 
attractor is 22.4±5.8 years (1 S.D.), for the CRI-HYM equilibrium is 31.1±6.7 years, 
and for the SEA equilibrium is 397.6 ± 308.2 years. 


6.6 Dependence of species interactions on neigh¬ 
borhood size 

The above results are based on the assumption that species interactions operate over a 
spatial scale of 600 cm 2 (i.e. the size of a single quadrat). Since the subtidal species in 
our model are sessile, however, they probably interact primarily with their neighbors. 
Thus species replacements are more likely to be effected by local neighborhood den¬ 
sities surrounding a patch than by the average population densities among quadrats 
(Pacala and Silander 1985, 1990). 

Here we examine the effect of spatial scale on species interactions by subdividing 
the quadrats into either 4 or 8 subsections. In effect, making the spatial scale of the 
quadrats smaller is equivalent to decreasing the neighborhood size (NS) around each 
patch. Subdividing the quadrats into 4 sections created 4 x 9 = 36 sub-quadrats, 
each with an area of 150 cm 2 . Subdividing the quadrats into 8 sections created 72 
sub-quadrats, each with an area of 75 cm 2 . Once the data was portioned in this 
manner, we followed the two step procedure outlined in section (6.4.3). 

Table 6.6 shows the best fitting models using a neighborhood size of 150 cm 2 
and 75 cm 2 . In all cases, the second best fitting model for each column of A x has 
a A AIC value greater than 4, suggesting that it has considerably less support than 
the best fitting model (Burnham and Anderson 1998). In general, decreasing the 
neighborhood size tends to reduce the number of species interactions in the model. 
For example, in the first column of A x (transitions from HYM), the Oy are a function 
of four interspecific interactions when NS = 150 cm 2 , but only three interspecific 
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Column of A 

Neighborhood Size 
150 cm 2 75 cm 2 

HYM 

m 4 

M 3 

MYX 

m 4 

m 2 

MYC 

m 4 

m 4 

SEA 

M 5 

Me 

CRI 

Mg 

M 5 

FIL 

m 4 

m 2 

FG 

m 4 

m 4 

BR 

M 6 

M 3 


Table 6.6: Best fitting models for each column of A obtained by dividing the quadrats 
into either 4 or 8 sub-quadrats. The sub-quadrats are equivalent to a neighborhood 
size of 150 cm 2 and 75 cm 2 , respectively. The symbol M h specifics the model with 
the lowest AIC value for the transition probabilities in each column of A, where the 
subscript h gives the optimal number of interspecific species interactions. 

interactions when NS = 75 cm 2 . The best fitting model for the 150 cm 2 neighborhood 
is M 44455446 and for the 75 cm 2 neighborhood is M 32455243 (Table 6.6). 

6.6.1 Equilibrium distributions 

The nonlinear Markov chain models M 44455446 and M 32455243 have multiple stable 
eq uili bri um states of attraction. Figure 6.3 shows the co-existing attractors at neigh¬ 
borhood sizes of 150 cm 2 and 75 cm 2 . In both cases the dominant attractor is the 
CRI-HYM equilibrium. Two other attractors coexist when NS = 150 cm 2 , an equi¬ 
librium distribution dominated by Filograna (FIL equilibrium) and an equilibrium 
distribution dominated by Mycale (MYC equilibrium). Both of these equilibria, how¬ 
ever, have much smaller basins of attraction than the CRI-HYM equilibrium. When 
NS = 75 cm 2 , the MYC equilibrium disappears completely, and the basin of the 
CRI-HYM attractor increases in size. 

To examine the effect of neighborhood size on the ability of the model to predict 
the species distribution in the subtidal community, we calculated the distance between 
the various CRI-HYM equilibria and the observed species distribution shown in figure 
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Figure 6.3: Coexisting attractors for the best fitting nonlinear Markov chain for 
neighborhood sizes of a NS = 150 cm 2 ), and b. NS = 75 cm 2 ). The percent in right 
hand corner of each histogram gives the proportion of times the model converged to 
each attractor starting from 10,000 random initial conditions. 














Model 

Distance ( D) 

Nonlinear 


NS = 600 cm 2 

0.68 

NS = 150 cm 2 

0.63 

NS = 75 cm 2 

0.44 

Linear 

0.19 


Table 6.7: Distance between the predicted and observed equilibrium distribution for 
nonlinear and linear Markov chains. Comparisons for the nonlinear Markov chains 
use the CRI-HYM equilibrium distribution. 


6.2. The distance D is measured as 

D = Y 1 \x i -6i\ (6.17) 

i=1 

where Xi is the predicted equilibrium frequency of species i and Oj is the mean fre¬ 
quency of species i. Table 6.7 shows values of D for neighborhood sizes of 600 cm 2 , 
150 cm 2 , and 75 cm 2 , and for the linear Markov chain M 0 (i.e. the are independent 
of species densities). While decreasing the neighborhood size in the nonlinear model 
decreases the distance between predicted and observed abundances, the linear model 
does a much better job at predicting community composition. 

6.6.2 Distribution of species interaction coefficients 

Figure 6.4 shows the distribution of species interaction coefficients for the best fitted 
models, estimated at different neighborhood sizes. In all cases the mean value is 
approximately 0, but the variance decreases with decreasing neighborhood size from 
105.6 (NS = 600 cm 2 .), to 44.2 (NS = 150 cm 2 ), to 19.52 (NS = 75 cm 2 ). Thus, as the 
neighborhood size decreases, the proportion of species interactions that are relatively 
weak increases. 
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Parameter Value 


Figure 6.4: Distribution of interaction coefficients for the best fitting nonlinear 
Markov chain estimated different size neighborhoods, a. NS = 600 cm 2 ). The total 
number of interaction coefficients is 371. b. NS = 150 cm 2 ). The total number of 
interaction coefficients is 308. c. NS = 75 cm 2 ). The total number of interaction 
coefficients is 266. 
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Equilibrium 

Neighborhood Size 

600 cm 2 150 cm 2 75 cm 2 

CRI-HYM 

FIL 

31.15 (6.73) 76.41 (26.78) 30.56 (26.28) 
— 9.61 (2.35) 13.19 (3.03) 


Table 6.8: Mean convergence time to the CRI-HYM and FIL equilibriums, based on 
10,000 simulations of the nonlinear Markov chain associated with each Neighborhood 
size. Values in parenthesis are 1 standard deviation (see text for details). 

6.6.3 Convergence times and patterns of succession 

The mean convergence times for the models M 44455446 (NS = 150 cm 2 ) and M 32455243 
(NS = 75 cm 2 ) to the CRI-HYM equilibrium and the FIL equilibrium are shown in 
Table 6.8. The mean convergence time to the CRI-HYM equilibrium is 3 to 8 times 
longer than the convergence time to the FIL equilibrium. There is no general relation 
between convergence times and neighborhood size. Convergence time for the CRI- 
HYM equilibrium is longest when NS = 150 cm 2 (intermediate spatial scale), while 
for the FIL equilibrium it is longest when NS = 75 cm 2 . 

Figure 6.5 shows the pattern of succession starting from an initial condition in 
which the abundance of bare rock is set to one, and abundances of all the species are 
set to zero. This represents succession following the total destruction of the existing 
co mmuni ty. When the neighborhood size is 600 cm 2 the community converges almost 
immediately to the periodic 2-cycle attractor (Fig. 6.5a.). When the neighborhood 
size is 150 cm 2 or 75 cm 2 , however, the community converges to the CRI-HYM equi¬ 
librium (Fig. 6.5b,c). Note that the time track of colonization differs between many of 
the species. For example, Filograna , the sea anemones, and the functional group (FG) 
reach an early peak in abundance within 2-4 years and then drop off fairly sharply 
before they approach their equilibrium. The sponges Hymedesmia and Myxilla , on 
the other hand, tend to increase monotonically towards their equilibrium abundance. 
Finally, Crisia increases monotonically when NS = 150 cm 2 , but shows an early peak 
in abundance followed by a slow decline towards equilibrium when NS = 75 cm 2 . 
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Species Frequency 



Time since disturbance (years) 


Figure 6.5: Simulation of state abundances over time, following complete destruction 
of the community (initial condition equals 100% bare rock), a. NS = 600 cm 2 ), b. 
NS = 150 cm 2 ), c. NS = 75 cm 2 ). 
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6.6.4 Rates of local succession 

A Markov chain describes two spatial scales—the dynamics of the community and 
the local dynamics of patches. At equilibrium the proportion of patches in each state 
is constant, however, individual patches change states continuously through time. 

Once the community converges to an equilibrium distribution the matrix A x re- • 
mains constant. From the equilibrium matrix A x we investigate rates of successional 
change at the spatial scale of a patch by calculating mean turnover rates and Smou- 
chowski recurrence times for each state. 


• The mean turnover rate describes the probability that a randomly selected patch 
chang es state between t and t + 1, and is given by 


^Wiil-au) (6.18) 

2=1 

where w { is the ith element of the dominant eigenvector, and (1 - an) is the 
probability that a patch in state i changes states in the time interval from t to 
t + 1 (Iosifescu 1980). 

• The Smouchowski recurrence time 0* of state i is the time elapsing between a 
patch leaving state i and then returning to it again. Its mean is given by 


0i 


1 ~ Wj 
ti?i(l O'ii) 


(6.19) 


(Iosifescu 1980). 

Table 6.9 shows mean turnover rates at the various equilibrium points. Patch 
turnover rates are much greater when the system has converged to the CRI-HYM 
equilibrium, than when it has converged the FIL or MYC equilibrium. In the case 
of the CRI-HYM and FIL equilibriums, the predicted turnover rate increases as the 
spatial scale of the density dependence decreases. 
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Equilibrium 

Density Spatial Scale 

300 cm 2 150 cm 2 75 cm 2 

CRI-HYM 

0.2700 

0.3112 

0.3655 

FIL 

— 

0.0031 

0.0786 

MYC 

— 

0.1673 

— 


Table 6.9: Mean turnover rate of patches when the subtidal community is at an 
equilibrium point. Turnover rates give the probability that a randomly chosen patch 
changes states in one time interval. The turnover rates in each column correspond 
to the best fit nonlinear Markov chain for a given spatial scale of density dependence 
(see text for details). 


Figure 6.6 shows Smouchowski recurrence times for each species at the CRI-HYM 
eq uili brium and the FIL equilibrium. In general, Qi decreases as the neighborhood 
size decreases. This is especially true in the case of the FIL equilibrium, where the 
recurrence time for many species decreases by as much as ten-fold when the NS 
decreases form 150 cm 2 to 75 cm 2 . 


6.7 Bifurcation Analysis 

When the neighborhood size of a patch is small, the nonlinear Markov chain always 
converges to a stable equilibrium. One of the signatures of nonlinear models, however, 
is their ability to produce periodic cycles, aperiodic oscillations and chaos. Shifts in 
parameter values can cause bifurcations from one of these behaviors to another (e.g., 
May 1974,1976, Hastings and Powell 1991; Pascual 1993, Costantino et al. 1995; 
Caswell and Neubert 1998). Here we look for bifurcations caused by changes in the 
strength of the interactions among species. We focus on the nonlinear Markov chain 
M 32455243 , parameterized using a neighborhood size of 75 cm 2 . 

The interaction coefficients for a given species can be changed in the following way. 
Let C ( N be a coefficient matrix whose (i,j) element is the interaction coefficient asso¬ 
ciated with species k in the function (x). The coefficient matrix for Hymedesmia 
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Figure 6.6: Smouchowski recurrence times for each patch state 0* at equilibrium. The 
different colored bars in the histograms represent recurrence times for different size 
neighborhoods, a. Smouchowski recurrence times for state i after the community has 
converged to the CRI-HYM equilibrium, b. Smouchowski recurrence times for state 
i after the community has converged to the FIL equilibrium. 





(k = 1) is shown below 

/ 

C (1) = 


V 

where the element is the value of the interaction coefficient associated with x x (i.e. 
the density of Hymedesmia ) in the function an(x), 4V is the value of the interaction 
coefficient associated with x\ in the function a 2 i(x), etc. The dashed symbols — 
indicate transition probabilities that are independent of the density of Hymedesmia. 
An example is the function a 32 (x) which represents the probability that Myxilla is 
replaced by Mycalle. 

To obtain a bifurcation diagram for species k we proportionally increase the in¬ 
teraction coefficients by a factor p 

C< fc) = C( fe > + (C« x p ) (6.21) 

and then insert the new coefficient values given in C£ fc) into the functions ay(x) to 
obtain a new transition matrix A x (fe,p) (the superscript (k,p) indicates that the inter¬ 
action coefficients associated with species k in the functions a tj have been increased 
by a factor of p). Bifurcation diagrams were obtained by starting the system at an 
equilibrium point and then increasing the value of p from 0 to 10 in increments of 
0.01. At each increment, we iterated the model 500 times to remove any transient 
dynamics, and then plot 50 iterates of the value of species k (i.e. x k ). 

Figure 6.7 shows a set of bifurcation diagrams starting from the CRI-HYM equi- 
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-2.87 

0.97 
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— 
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librium, in which the value of Xf. is plotted as a function of the proportional change 
(p) in the interaction coefficients of species k. Increases in the interaction strength of 
Crisia and Mycale produce the most dramatic bifurcation patterns, including periodic 
cycles, phase locking and chaotic dynamics. Increases in the interaction strength of 
Myxilla results in a bifurcation from a stable equilibrium to a 2-cycle at p = 3.67, 
while increases in the interaction strength of Filograna and the sea anemones produce 
a sharp transition in their frequencies (from less than 0.05 to 1) at p equals 3.43 and 
0.78, respectively. For the remaining species, increasing p results in a smooth change 
in eq uili bri um abundance, but no bifurcations. 

The most dramatic bifurcation patterns are associated with Crisia (Figure 6.8). 

If we start the system in the basin of attraction of the CRI-HYM equilibrium (Fig. 
6 .8a), the attractor is a fixed point for p < 2.67. At p — 2.67 a Hopf bifurcation occurs 
in which the fixed point is replaced by an invariant circle. At p ~ 3.51, the invariant 
circle is replaced by a 3-cycle. Increasing p further results in period doublings and 
eventually leads to chaos. In the chaotic region of parameter space, species densities 
are aperiodic in time and show sensitive dependence to initial conditions (Fig. 6.9). 
When we start the system in the basin of attraction of the FIL equilibrium, the 
attractor remains a stable fixed point until p « 7.09 (Fig. 6.8b). At p > 7.09 the 
fixed point becomes unstable and the system falls into the basin of attraction of the 
chaotic attractor. 

Increasing the interaction strength of the remaining species has no effect when we 
start the system in the basin of attraction of the FIL equilibrium. Thus even though 
the basin of attraction of the FIL equilibrium is comparatively small (relative to the 
CRI-HYM equilibrium), the FIL equilibrium is more stable to large perturbations in 
interaction strengths. 

Because many interaction coefficients are positive, increasing the interaction strength 
of one species increases the relative abundance of other species in the community. Fig¬ 
ure 6.10 shows how the abundance of Hymedesmia , Mycale , and the FG species vary 
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HYM 


MYX 



Proportional increase (p) in the interaction 
coefficients of state k 


Figure 6.7: Bifurcation diagrams plotting the frequency of each species k (k = 
1 ,..., 8) as a function of the proportional increase (p) in the interaction coefficients 
of species k (NS = 75 cm 2 ). Each bifurcation diagram was generated by starting from 
an initial condition in the basin of attraction of the CRI-HYM equilibri um . 
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Figure 6.8: Bifurcation diagrams plotting the frequency of Crisia on the attractor 
as a function of the proportional increase (p) in the interaction coefficients of Crisia 
(NS = 75 cm 2 ), a. Bifurcation diagram starting on the CRI-HYM attractor, b. 
Bifurcation diagram starting on the FIL attractor. 
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Figure 6.9: Sensitivity to initial conditions. Time series plot of Crisia showing that 
temporal solutions diverge for small initial differences. Two trajectories of Crisia are 
shown in which the initial frequency of Crisia differs by only 0.001 (p = 8). 
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as a function of the proportional increase in the interaction strength of Hymedesmia. 
Increasing p from 0 to 1 increases the abundance of Hymedesmia. At p ~ 1, however, 
there is an abrupt decline in Hymedesmia accompanied by a sharp increase in the 
FG species. When p > 1, Hymedesmia abundance is a decreasing function of p while 
Mycale abundance is a monotonically increasing function of p. Thus, increasing the 
interaction strength of Hymedesmia enables Mycale and the FG species to increase 
their abundance at the expense of Hymedesmia. 


6.8 Discussion 

There have been very few empirical studies that have tried to quantify interaction 
strengths among assemblages of species. Most studies have relied on estimating per 
capita interaction strengths through single species manipulation experiments (e.g., 
Fowler 1981, Menge and Farrell 1989). The problem with this approach is that there 
may be a large number of interactions and that the prevalence of indirect effects 
in most communities (the effects of one species on another via its effect on a third 
species) means that pair-wise estimates of interaction strengths are not additive. 
Thus, extrapolating results from single species manipulation experiments to whole 
communities are likely to be unsuccessful (Wilbur 1972; Neil 1974; Wilbur and Fauth 
1990; Wooten 1993). 

The methods presented here allowed us to estimate density-dependent interac¬ 
tions in an unmani pulated community using spatial time series data. By focusing on 
patch tr ansi tion probabilities, instead of individual species, we were able to simul¬ 
taneously estimate interaction coefficients for all species states in our model using 
maximum likelihood methods. Any indirect effects are automatically incorporated 
into the model because interaction coefficients are estimated by taking into account 
the dynamics of the community as a whole. This is not true for species removal 
experiments. 
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Proportional increase (p) in the interaction 
coefficients of HYM 


Figure 6.10: Bifurcation diagrams for the abundance of the states HYM, MYC, FG 
and BR as a function of the proportional increase (p) in the interaction coefficients 
of Hymedesmia. 
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Another difference between our methods and species removal experiments is that 
removal experiments tend to focus on interesting or significant interactions, which 
means they are usually selected because they axe strong (Lawton 1992). Our analy¬ 
sis shows not only that the best fitting models contains a large number of species 
interactions but also that most of these interactions are relatively weak. This re¬ 
sult is consistent with several empirical studies of plant communities showing they 
are dominated by a large number of competitively weak interactions (e.g., Schoener 
1983, Mahdi et al. 1989, Polis 1991). 

Parameter estimation: The two-step method outlined in section 6.4.3 allowed 
us to identify the best fitting nonlinear Markov models, within an extremely large 
parameter space, in a relatively short period of time (2 to 4 days depending on the 
size of the neighborhood). The reason for this is two fold. First, the binary logistic 
regression method greatly reduces the number of species combinations one needs to 
test to find the best fit for a given model. For example, to find the best M x model for 
column j, one needs to identify the species whose density has the greatest effect on 
the replacement of species j by species i, i = 1 ,..., 5 (i.e. the highest ranked species 
for each ay(x)). Since the number of possible choices per transition probability is 
5—1, the n umb er of possible species combinations we need to test to find the highest 
ranked species for all the transition probabilities in column j is 5(5 — 1). In contrast, 
if we used polychotomous logistic regression to perform the same task, the number of 
possible combinations we would need to test is (5 — l)*- 5-1 ). The problem, of course, 
only gets worse as the number of species interactions increases. 

The second reason the two-step method shortens the search time is that maxi¬ 
mization of Eq. 6.11 (binary likelihood function) is much faster than Eq. 6.6 (poly¬ 
chotomous likelihood function). Thus, the binary regression methods enable us to 
quickly identify the relative ranking of the species, according to their effect on the 
replacement of species j by species i, for each entry in A x . Once the species rankings 
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axe determined, however, the interaction coefficients associated with a given model 
M h axe found fairly quickly using the polychotomous regression method (since we 
only need to maximize Eq. 6.6 once for a given set of species interactions). 

Community dynamics: Simulation of the nonlinear Markov chain suggested that 
the subtidal community is capable of periodic orbits and multiple eq uili bria. This 
suggests that large perturbations in species abundance (caused by random climate 
change or biological disturbances) may result in ecological shifts from one eq uili bri um 
region to another. This pattern of multiple equilibrium states has been observed 
in several empirical studies of natural communities (Holling 1973, 1992; Sutherland 
1974; Levin 1976; May 1977; Hughes 1994). Sutherland (1974) and Connell and 
Slatyer (1977) have extensively reviewed work on the marine rocky intertidal, on 
coral reefs, on freshwater lakes, and on terrestrial plant communities, and concluded 
that community structure is indelibly tied to historical events and therefore multiple 
equilibrium states are an undeniable reality of natural systems. 

Our model of the subtidal community corroborates these conclusions and suggests 
that large-scale phase shifts in community composition are possible. If we assume 
that the current state of subtidal communities in the Gulf of Maine falls within the 
basin of attraction of the CRI-HYM equilibrium, then any change in the physical 
environment, which pushes the system out of that basin, will drastically change the 
structure of the community. All of the other equilibrium states in our model predict 
complete domination of the substrate by a single species (regardless of the density 
spatial scale used to parameterize the model). Thus perturbations that shift the 
community towards a new equilibrium will drastically reduce the diversity of the 
community. Huges (1994) has documented such an event in coral reef co mmuni ties off 
Jamaica, where overfishing, disease, and hurricane damage have combined to destroy 
most corals, resulting in dramatic phase shift to a low diversity system dominated by 
fleshy macroalgae. 
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Distribution of interaction coefficients: We found that interaction strengths in 
the subtidal community were symmetrically distributed around zero, with interaction 
coefficients equally as likely to be positive as negative. While negative coefficients 
probably represent competition between species for space, positive coefficients repre¬ 
sent not only mutualistic interactions but also indirect species effects (e.g., Callaway 
1997, Hodge et al. 1999, Levine 1999). If the presence of species X , inhibits species Y 
which has a negative effect on species Z, then the interaction coefficient corresponding 
to species X in the function a Zj (x.) is likely be positive. Paine (1992) found similar 
indirect positive interactions between intertidal invertebrates (chitons and limpets) 
and two species of brown algae ( Hedophyllum sessile and Alaria ssp.). In this case 
the invertebrates prey on a crustose coralline algae, which in turn inhibits the recruit¬ 
ment of H. sessile and Alaria. Thus, the actions of the invertebrate predators have 
a positive effect on the brown algae because it increases the recruitment rate of H. 
sessile and Alaria onto the substrate surface. 

The effects of neighborhood size: The number of species interactions and the 
distribution of interaction strengths in our model is dependent on neighborhood size 
we chose to parameterize the model. When density effects are measured over smaller 
spatial scales the number of interactions decreases and interaction strengths became 
weaker (Fig. 6.4). The pattern of species interactions specified by our model is 
similar to several recent empirical community studies which show that most species 
have weak effects on the abundance of other species, while only a few have strong 
effects (Paine 1992, Fagan and Hurd 1994, Raffaelli and Hall 1996; Wootton 1997). 

The neighborhood size used to estimate species densities also affected the behavior 
of the nonlinear Markov chain. When we parameterize the model using an neighbor¬ 
hood size of 600 cm 2 , the dominant behavior of the system is a two-cycle, in which the 
co mmuni ty is dominated by the sponge Mycale. As we decreased the neighborhood 
size, however, the best-fitting model predicts a more stable and diverse community 
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structure (CRI-HYM equilibrium). This may be due to the increase in the number of 
weak interactions between species. Weak species links have been shown to dampen 
oscillations in theoretical food web models, thus reducing the probability that species 
go extinct (McCann et al. 1998). 

Finally, neighborhood size also plays an important role in measurements of local 
patch dynamics. Mean turnover rates and species recurrence times both increase 
as the neighborhood size decreases. This result points out one of the fundamental 
dilemmas in ecology; identifying the characteristic scale at which local interactions 
operate to produce observed community patterns at larger temporal and spatial scales 
(Levin 1992; Pascual and Levin 1999). 

Nonlinear Dynamics: The bifurcation analysis allowed us to examine what would 
happen to the community if the intensity of interactions associated with one species 
were increased. Increases in the strength of the interaction of the bryozoan Crisia 
(and to a lesser extent Mycale) with the rest of the community lead to periodic oscil¬ 
lations and chaotic dynamics. While chaos has been found in numerous theoretical 
community models (e.g., Gilpin 1979, Hastings and Powell 1991, Kot et al 1992, 
Caswell and Neubert 1998), this is the first time it has been reported in a Markov 
chain model for community competition. 

The bifurcation analysis also shows that increasing the interaction strength of a 
species does not guarantee its abundance will increase. For example, increases in the 
strength of the interaction of Hymedesmia with the rest of the co mmuni ty actually 
lead to a decrease in its abundance (Fig. 6.10). This occurred because Hymedesmia 
has a strong negative effect on the abundances of Crisia, Myxilla, Filograna, and the 
sea anemones, which in turn have a negative effect on the abundance of Mycale. If the 
interactions of Hymedesmia with other species become too strong, then the species 
inhibiting Mycale are almost completely wiped out. This allows Mycale to increase 
its abundance in the community and actually inhibit Hymedesmia. 
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6.9 Conclusion 


Most applications of Markov chain models to marine systems have been limited to 
linear models, in which transition probabilities are time-invariant. Patch transitions, 
however, certainly depend on species abundances. To assess the importance of such 
factors in a rocky subtidal community, we constructed a nonlinear Markov chain by 
making patch tr ans ition probabilities a function of species densities. We introduced a 
two-step method, based on logistic regression analyses for binary and polychotomous 
variables, that allowed us to estimate the density dependence of the elements A x . 
Parameters for these functions were estimated using maximum likelihood analysis 
and the best models chosen using AIC methods. 

By applying this methodology to the subtidal data set, we showed that the commu¬ 
nity possess multiple equilibrium states, that the distribution of species interactions 
are symmetrically distributed around zero (and mostly weak), and that the behavior 
of the model depends on the spatial scale (neighborhood size) at which the parame¬ 
ters are estimated. These results illustrate the complexities involved in predicting 
the behavior of any species assemblage. First, systems with multiple equilibria can 
change radically given a large perturbation (May 1977) and thus future states of the 
community are highly dependent on random environmental events. Second, choosing 
the correct spatial scale to parameterize the model is not a simple task, and nonlinear 
models parameterized using different density spatial scales can give entirely different 
results (Pascual and Levine 1999). 

Our approach diff ers from traditional methods of measuring interaction strengths 
because we are concerned with their effect on patch transitions as opposed to species 
abundances. The advantage of our methodology is that it allows us to simultaneously 
estimate a set of interaction coefficients based solely on time series data. Because of 
the no nlin earities inherent in ecological systems, the effect of changes in interaction 
strengths and/or species densities can only be understood if all possible combina¬ 
tions of interactions are estimated at the same time (Bender et al. 1984). Pair-wise 
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estimations of interaction strengths (via species removal experiments) provide little 
information about how species will interact when they co-occur within an assemblage 
of organisms, and what effect those interactions will have on the dynamics of the 
community (Wootton 1993). 
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Appendix 


The tables below show the rank ordering of species interaction strengths for each of 
the probability functions a^- (x) pertaining to the best-fit model (Neighborhood Size 
= 600 cm 2 ). Each row in the tables identifies the species whose densities affect the 
transition from state j to state i. The order of the species in the rows gives the species 
ranking, from strongest to weakest, according to their affect on the probability that a 
patch in state j changes to state i in one time interval (see section 6.4.3 for details). 
The species in the columns designated as ” Out” are not included in the probability 
functions (x). 
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Chapter 7 


Conclusion and Future Directions 


In this thesis I have presented a series of spatial models describing the dynamics of 
metapopulations and a rocky subtidal community. These models are motivated by 
recent theoretical and empirical evidence that the distribution of organisms in space 
can have important consequences for the function and structure of terrestrial and 
aquatic systems (Levin and Paine 1974; Steel 1978; Pickett and White 1985; Caswell 
and Cohen 1991a,b; Weins et. al. 1993; Wu and Levin 1994; Wu and Loucks 1995; 
Hanski 1999). 

While the questions addressed in each of the chapters are diverse, the underlying 
assumption of the models is that an assemblage of organisms can be represented as 
a mosaic landscape of discrete patches. Patches can be in one of N possible states, 
where the state of a patch is defined by the organism or organisms occupying it. In 
the metapopulation models (Chapters 2 and 3) patches represent habitat fragments 
capable of supporting a local population. In subtidal community models (Chapters 
4, 5, and 6) a patch represents the spatial location of an invertebrate species on the 
rock wall substrate. In either case, however, the models are formulated in terms of a 
transition matrix whose entry gives the probability that a patch in state j changes 
to state i (i = 1,..., N ) in a single time step. 
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7.1 New insights into the effects of habitat de¬ 


struction 

In the past ten to fifteen years, metapopulation models have become an important 
tool for studying the effects of habitat destruction on metapopulation persistence. 
While modeling methods have varied, the general thrust of these studies has been to 
predict how habitat destruction will affect occupied patch frequencies and to quantify 
the amount of habitat destruction the population can tolerate before it goes extinct 
(the so called extinction threshold; Tilman et al. 1994). Today the identification of 
extinction thresholds for endangered species has become an important conservation 
strategy. 

Metapopulation models have typically been predicated on three assumptions 

1. The landscape is composed of an infinite number of habitat patches. 

2. Propagules are capable of dispersing throughout the landscape; i.e. colonization 
of an empty patch is not affected by distance. 

3. Suitable patches are randomly distributed in space. 


Extinction thresholds on fractal landscapes: Chapter 2 looks at the effect of 
relaxing assumptions 2 and 3, by exploring the effects of habitat destruction in a 
spatially explicit cellular automaton (CA) model. In the CA, the dispersal range 
of propagules is limited to a 3-patch radius and the habitat destruction pattern is 
explicitly defined using fractal maps. The general findings of this study are illustrated 
in Figure 7.1. When dispersal is local, the equilibrium frequency of occupied patches 
(p) decreases and the extinction threshold (h c ) increases as the fractal dimension 
( D) of the habitat destruction pattern increases. The take home message is that in 
large landscapes (consisting of 100’s if not 1000’s of patches), the spatial structure of 
suitable habitat is at least as important as the amount of suitable habitat. 
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Extinction Threshold 


Figure 7.1: Generalized plots of the equilibrium frequency ft and the extinction thresh¬ 
old h c as a function of the fractal dimension of the landscape. The arrows on the right 
and left of the graphs show that decreasing the fractal dimension of the landscape 
is equivalent to increasing the connectivity of suitable patches and decreasing the 
amount to edge between suitable and unsuitable habitat. 
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In today’s modern world, habitat loss is inevitable. While the situation is grim for 
many endangered species, efforts should at least be made to ensure that habitat de¬ 
struction patterns have the lowest possible fractal dimension. The fractal dimension 
of a landscape can be easily measured (Krummel et. al. 1987; Sole and Manrubia 
1995; Milne 1988) and conservation efforts to minimize D will maximize connectivity 
in the remaining suitable territory while minimizing edge effects (Fig. 7.1). Lowering 
the fractal dimension of habitat destruction patterns may at least give some endan¬ 
gered species a fighting chance. Perhaps this is the best we can hope for in a globalize 
economy bent on extracting resources at an ever increasing pace. 

Chain-Binomial Metapopulation Model: Chapter 3 looks at the effects of re¬ 
laxing assumption 1, by exploring the effects of habitat destruction in finite patch 
networks. The chain-binomial metapopulation (CBM) model develop to describe this 
situation is applicable to marine populations that produce planktonic larvae capable 
of long distance dispersal. In finite landscapes there is no such thing as a posi¬ 
tive equilibrium frequency or extinction threshold, because with probability 1 the 
metapopulation always goes extinct. The important problem for such systems is to 
determine how habitat destruction affects the extinction time, r, of the metapopula¬ 
tion. The CBM model reveals a dangerous scenario regarding attempts to predict the 
effects of habitat destruction — f declines greater than exponentially as the number 
of suitable patches, S, declines. Thus, a small amount of habitat destruction has 
the potential to drastically reduce t from effectively infinite values (on an ecological 
time scale) to very small values, leaving a once seemingly healthy population on the 
verge of extinction. The sensitivity of f to changes in S makes the task of identifying 
critical habitat thresholds in finite landscapes a near impossibility. This result has 
especially grave consequences for endangered species, who by definition are confined 
to small patch networks. 
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7.2 Markov chain models of sessile communities 


The structure of marine communities is an emergent property of species interactions 
occurring over relatively small spatial scales and physical processes occurring over a 
range of spatial scales. Determining how these processes interact to produce complex 
community structures, however, remains a fundamental problem in marine ecology. 

In the second half of this thesis, I have shown that Markov chains have the po¬ 
tential to vastly improve our ability to identify key ecological factors (succession, 
disturbance, competition, environmental variability, etc.) and key species that are 
important for the structure and function of marine sessile communities. The ap¬ 
proach I have taken is to start with simple Markov chain models, which exclude 
much of the detail of the system, before proceeding to more complex models. Even 
the most complicated of them, however, is not intended as a detailed quantitative 
description of all the factors affecting rocky subtidal communities. My philosophy 
is that comparing related models that differ in biological details is a more powerful 
approach for understanding community dynamics than analyzing any one model. 

Linear Markov chain models: In chapter 4 I began the analyses of the rocky 
subtidal co mmuni ty by constructing a linear time-invariant Markov chain. This is the 
simplest model that incorporates the basics of patch transitions, species replacement, 
and disturbance. While the linear Markov chain ignores the effects of environmental 
variability and density dependence, it accurately predicted species distribution in the 
Gulf of Maine, provided ecological information for categorizing species into functional 
groups, and revealed that 38% of the substrate was occupied by a different species 
each year. 

One of the biggest advantages of linear Markov chains is that sensitivity analysis 
can be developed to determine how changes in model parameters affect predicted com¬ 
munity patterns. Sensitivity analysis has already become an essential part of demo¬ 
graphic analysis (Tuljapurkar and Caswell 1997, Caswell 2000), and has the potential 
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to become equally valuable in community modeling. Here I used sensitivity analysis 
to characterize how changes in elements of the transition matrix affect predictions 
of species abundance and community diversity. The results suggest that diversity 
is much more sensitive to changes in colonization probabilities than to species re¬ 
placement probabilities. This finding is consistent with observations by Witman and 
Sebens (1990), who found that competition plays a relatively minor role in structuring 
subtidal communities. 

Linear Markov chain models are also ideal for quantifying the effects of species 
removal on diversity and community resilience. Analysis of the subtidal transition ma¬ 
trix showed that the removal of Hymedsmia spp. 1 or Myxilla fimbriata produces the 
largest chang e in diversity, however, the overall change is relatively small. This find¬ 
ing indicates that there axe no real keystone species in the subtidal community. The 
removal of bare rock or Crisia ebumea, however, produced relatively large changes 
in the rate of converge to equilibrium, suggesting that these states are important for 
community resilience. The removal analyses developed in chapter 4 are Markov chain 
analogs to species removal experiments in natural systems. They provide a means of 
classifying the relative importance of individual species to the structure and stability 
of the community when experimental manipulations are not possible. 

Environmental variability: The next step in the complexity hierarchy (chapter 5) 
was to test for significant effects of temporal and spatial variation on patch transi¬ 
tion probabilities. Because this is a 2-factorial test, I used log-linear analysis, which 
showed that transition probabilities varied significantly between quadrat locations 
and between years in the subtidal zone. While some authors (e.g. Li 1996, Childress 
et al. 1998) have argued that such variability must be taken into account for Markov 
rhains to adequately predict equilibrium community structure, I found that including 
this variability in the subtidal model had little effect on predicted species abundances. 
On the other hand, simulation of the time-varying Markov chain shows that temporal 
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variation in successional processes (i.e. transition probabilities) has important conse¬ 
quences for transient dynamics. In the subtidal community the time-varying model 
predicts frequent shifts in composition between a Hymedesmia dominated co mmu nity 
and a Crisia dominated community. This pattern is consistent with observations in 
the Gulf of Maine over the past ten years (Witman 1996). 

Chapter 5 also introduces a method for comparing probability matrices (Cohen 
et al. 1998) and uses it to characterize regional variability in the distribution of 
species in space and time. While this methodology is somewhat of a digression from 
the Markov chain approach, it provides a simple but effective way of characterizing 
the spatial and temporal stability of individual species relative to other organisms in 
the community. In the subtidal community, Cohen’s method predicts there is more 
variability between quadrats in the temporal distribution of species, than there is over 
time in the spatial distribution of species. This result is consistent with the log-linear 
analysis, which shows a greater effect of space on transition probabilities than time. 

Density dependence Chapter 6 increases the complexity in the subtidal models fur¬ 
ther by incorporating the effects of density dependence on transition probabilities. 
The most important contribution of this work is the development of maximum likeli¬ 
hood methods to estimate the effects of species densities on transition probabilities. 
The likelihood methods will allow ecologist to estimate density-dependent interac¬ 
tions in an unmanipulated community from spatial time series data. By focusing 
on patch transition probabilities, instead of individual species, the strength of inter¬ 
actions among all species can be estimated simultaneously. Any indirect effects are 
automatically incorporated into the model because species interactions are estimated 
by taking into account the dynamics of the community as a whole. This is not true for 
species removal experiments, which have traditionally been used to measure species 
interaction strengths in community assemblages (Paine 1992). 
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7.3 Future research directions 


Below I briefly describe some future research directions related to the modeling meth¬ 
ods presented in this thesis. 

7.3.1 Metapopulations models: Habitat restoration and Ma¬ 
rine reserves 

Habitat restoration: To study the effects of habitat destruction I assumed that 
once habitat was destroyed it remained destroyed forever. In many marine systems, 
however, there is a continuous turnover of suitable habitat, in which the destruction 
of habitat in one region is balanced by the restoration of habitat in another. A classic 
example is hydrothermal vents (Mullineaux et al. 1991), in which subpopulations 
among vent sites can be described as a metapopulation living on a geographical ex¬ 
panse of ephemeral patches. Extinction of a subpopulation results when the vent site 
it occupies becomes inactive (a process equivalent to habitat destruction). Propag- 
ules produced by many vent species, however, can disperse over long distances (Lutz 
1988) and are capable of colonizing newly developing vent sites (a process equivalent 
to habitat restoration). 

To study the characteristic of such systems I will include habitat restoration in 
the metapopulation models presented in chapter 2. Figure 7.2 shows an idealize tran¬ 
sition diagram for a metapopulation model with habitat restoration. The goal of this 
research is to analyze the system of equations specific by the transition diagram in 
order to characterize how the rate of destruction and recovery affect persistence and 
extinction times. I will then translate this model into a spatially explicit cellular au¬ 
tomaton model to study how the spatial location of habitat sites affects the dynamics 
of the system. 
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Suitable habitat 



Figure 7.2: Generalized transition diagram for a metapopulation model with habitat 
destruction and habitat restoration 
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Non-Reserve 


Reserve 




Figure 7.3: Transition diagram for a hyper-population model of a marine reserve. 
The proportion of the area that falls within the reserve equals q, where 0 < q < 1. 
Xi and X 3 represent unoccupied patches in the non-reserve and reserve, respectively. 
X 2 and X A represent occupied patches in the non-reserve and reserve. Q (i = 1,2) is 
the probability that an unoccupied patch is colonized. (5* (i = 1,2) is the probability 
that a local population goes extinct. See text for details. 

Marine reserve design: Setting aside reserve areas for protection from fishing 
is an important option being considered in the management of fisheries (Clark 1996; 
Fogarty 1999). To study how marine reserves could help conserve fish populations and 
benefit fisheries I have developed a hyper-population model describing the dynamics 
of a pair of metapopulations located in two distinct regions. Figure 7.3 shows a 
transition diagram for the hyper-population model, in which the distinct regions 
correspond to reserve and non-reserve areas in a fishery. Fishing is only allowed in 
the non-reserve areas, and the proportion of the fishery that is set aside as non-reserve 
is 1 — q (0 < q < 1). 

Metapopulations in the reserve and non-reserve regions are connected to each 
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other via the colonization functions C\(x 2 , x 4 ) and C 2 (£ 2 , x 4 ), given by the equations 

n ( \ -i fb(x 2 (l-a 2 i) + x 4 a 12 )\ 

C i(x 2 , x 4 ) = 1 - exp I- - -J 

n , \ , /6(x 2 a 2 i + &2Z4 (1 - ai 2 ))\ , . 

C 2 (a; 2 , 2 : 4 ) = 1 - exp I-—--I (7.1) 

where b is the reproductive parameter (see Eq. 2.10), Xi {1 — 2,4) is the proportion 
of occupied patches in the fishery, a i2 is the proportion of propagules produced in 
the reserve that end up in the non-reserve (0 < Qi 2 < q), and q 2 i is the proportion of 
propagules produced in the non-reserve that end up in the reserve (0 < a 2 i < 1 — q). 
Thus, the parameters o;i 2 and a 2 \ define the level of mixing of propagules between 
reserve and non-reserve areas. 

Fishing mortality in the non-reserve areas increases the disturbance rate of occu¬ 
pied patches and the total fishing pressure is assumed to be constant (i.e. fisherman 
don’t stop fishing just because some of the fishing grounds are set aside as reserves). 
This means that the fishing pressure per unit area increases as q increases. This effect 
can be modeled as 

( f 

81 = 1 — exp — [ k + - - 

V 1 ~Q 

S 2 = 1 — exp (k) (7.2) 

where k is a disturbance parameter and / is a measure of the total fishing pressure. 

The marine reserve metapopulation model is applicable to benthic species, such 
as scallops and oysters, as well as coral many reef fish (Man et al. 1995). I propose 
to analysis this model to determine the optimal reserve size for maximizing the ex¬ 
ploitable stock and the sustainable yield of the fishery. I will study how this optimum 
varies as a function b, k, /, an and a 2 i. I will also translate the reserve-model into a 
spatially explicit cellular automaton to look at the effect of the spatial arrangement 
of reserve areas on stock size and maximum yield. 
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7.3.2 Spatially explicit models of benthic communities 


Spatial-Dependent Markov Chain Models: In this project I will develop a spa¬ 
tially explicit model of the subtidal community by identifying patch transition proba¬ 
bilities that are dependent on neighboring species densities. This approach is similar 
to the methods I used in chapter 6 except that the likelihood functions must estimate 
species interaction effects on a patch-by-patch basis. Thus, each cell might have its 
own vector of independent variables, and the parameters would be estimated over the 
entire set of observed transitions. This approach will be computationally intensive, 
but it is worth exploring because it may give different results from the previous meth¬ 
ods, which only approximated spatial-dependent effects by dividing the quadrats into 
smaller subunits. 

I will use the information from this analysis to build a cellular automaton (CA) 
model in which patches transitions are a function of the states of neighboring patches. 
There is wide latitude in the form of the neighborhood and the nature of the transition 
functions. The difference between a CA model and the nonlinear Markov chain model 
is that the explicit spatial arrangement will affect the outcome. I will use the CA 
model to explore these effects on coexistence and spatial pattern, and explore the 
possibility of long-term phase shifts in community composition. 

A Mechanistic Cellular Automaton Model: State transitions are actually de¬ 
termined by mechanisms that involve growth, recruitment, disturbance, and com¬ 
petition for space. Here I propose to develop a cellular automaton model based on 
measurements of these mechanisms. My goal is to produce a matrix that will give the 
probabilities of transitions among states as a function of the state of the neighboring 
patches. 

The procedure will be as follows. 
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1 . First, I will estimate the probabilities of each of the basic interactions. My goal 
is to obtain estimates that give, as nearly as possible, the probabilities of each 
process in the absence of the others. 

(a) I will estimate the probability that an empty patch is colonized by each 
species from a set of species removal experiments (Witman, unpublished 
data). These experiments were performed by clearing 30 x 20 cm 2 areas 
of all organisms and following recolonization over a two year period. 

(b) I will estimate probabilities of growth into empty patches from the pho¬ 
tographic records of quadrats over time. This will require some ingenu¬ 
ity, because growth rates (which are directly measurable) depend on the 
amount of available space and the species present in neighboring patches. 

(c) I will estimate probabilities of competitive overgrowth by examining loca¬ 
tions where two species occupy adjacent patches. 

(d) I will estimate probabilities of disturbance by monitoring the appearance 
of empty patches. I will use the method of Caswell and Etter (1993) to 
model disturbances of different sizes. 

2. The different processes by which a patch can change state are a set of competing 
risks (e.g., Chiang 1966, David and Moeschberger 1978). Twill use competing 
risk theory to calculate the probabilities of changing state in the presence of 
competing risks. These probabilities will form the cellular automaton transition 
matrix. 

I will study the model by simulation, in order to determine how changes in recruit¬ 
ment, growth, and competitive exclusion rates affect community structure. I will also 
investigate the impact of the frequency and size distribution of disturbance events on 
model behavior. Finally, I will compare predictions of the mechanistic CA with the 
non-mechanistic models in order to determine which method better characterizes the 
dynamics of the rocky subtidal community. 
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